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Since  altimeter  measurements  are  directly  affected  by  the  surface  (geocentric) 
tide,  an  exact  representation  of  the  latter  is  important.  An  improvement  of  the 
altimeter  model  with  tidal  effects  included  is  achieved  by  an  adaptation  of 
Schwiderski 's  formula  giving  the  ocean  bottom  deformation  due  to  ocean  tidal  load-J 
ing.  The  resulting  model  for  the  bottom  tide  and,  especially,  for  the  surface 
tide  can  be  used  in  conjunction  with  all  of  the  tidal  constituents. 

Due  to  an  imperfect  description  of  the  earth's  gravity  field  by  a  truncated 
spherical-harmonic  series  in  a  first-phase  adjustment  process,  thei/neglected 
effect  is  superimposed  on  the  altimetric  noise  and  acts  as  if  increasing  its 
variance  and  introducing  covariances  in  a  global  sense.  This  added  effect  can 
be  evaluated  using  the  geoidal  covariance  function.  The  total' effect  (including 
the  noise  proper),  when  treated  as  the  first-order  autoregressive  process,  yields 
a  highly  patterned  variance-covariance  matrix  whose  inverse  is  a  tri-diagonal 
matrix.  The  latter  is  used  here  as  the  weight  matrix  of  observations  in  a  least- 
squares  adjustment,  offering  almost  the  same  computer  economies  as  a  diagonal 
weight  matrix  when  compared  to  the  rigorous,  fully  populated  weight  matrix.  How¬ 
ever,  the  errors  associated  with  a  diagonal  weight  matrix  are  greatly  reduced  in 
this  approach. 

^The  point-mass  adjustment  model,  based  on  the  Residuals  from  the  spherical- 
harmonic  adjustment  of  satellite  altimetry,  -ha^&een  recently  modified  to  allow 
for  an  efficient,  large-scale  resolution  of  the  geoidal  detail.^  That  task  has 
been  accomplished  with  the  aid  of  an  algorithm  resulting  in  a  banded  system  of 
normal  equations.  However,  the  point-mass  model  has  not  included  any  tidal 
effects.  In  the  present  analysis  this  model  is  extended,  in  the  sense  that  cer¬ 
tain  tidal  parameters,  associated  with  a  given  ocean  basin,  are  adopted  as  adjust-| 
able  quantities.  In  a  least-squares  adjustment  encompassing  such  a  basin  the  new 
set  of  parameters  is  common  to  every  observation  and  leads  thus  to  a  full  border 
portion  of  the  matrix  of  normal  equations.  This  results  in  a  banded-bordered 
system  of  normal  equations  which  can  be  resolved  almost  as  efficiently  (due  to  a 
relatively  small  number  of  tidal  parameters)  as  the  above  banded  system. 

■Us  an  alternative  to  the  point-mass  adjustment  (without  tidal  parameters),  an 
approach  is  presented  which  is  an  adaptation  of  the  collocation  theory.^First, 
the  standard  prediction  formulas  corresponding  to  the  least-squares  collWation 
with  noise  are  recapitulated,  yielding  the  "unsmoothed"  geophysical  quantifies 
based  on  the  unaltered  signal  and  noise  at  every  observation  point.  Subsequently, 
these  formulas  are  modified  to  yield  "smoothed-out"  geophysical  quantities.  Here 
a  part  of  the  signal,  representing  the  (fine)  gravity  field  variations  beyond  the 
desired  smoothing  level,  is  pushed  into  the  realm  of  "noise".  The  standard  col¬ 
location  approach  is  unchanged  except  that  during  the  formation  of  the  first 
matrix  in  the  prediction  formula,  the  covariance  (or  cross-covariance)  function 
is  no  longer  utilized  in  conjunction  with  the  expansion  degree  but  only  in 
conjunction  with  the  degree  corresponding  to  the  required  smoothness.  Finally, 
the  (pure)  densification  of  the  predicted  geophysical  quantities  in  view  of  con¬ 
tour  maps  can  be  carried  out  efficiently  via  the  errorless  collocation.  Con¬ 
sistent  with  this  procedure,  the  residuals  at  observation  points  computed  with 
respect  to  the  geoidal  contour  map  can  also  be  obtained  using  the  errorless 
collocation. 
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1.  INTRODUCTION 


This  report  pursues  the  development  of  satellite  altimeter  adjustment 
model  based  on  the  short-arc  algorithm.  In  the  first-phase  gravity  field 
adjustment  the  data  consist  of  altimeter  observations,  the  smoothed  surface 
approximating  the  geoid  is  described  by  a  truncated  set  of  spherical- 
harmonic  potential  coefficients,  and  each  satellite  arc  is  described  by  six 
state  vector  components.  In  addition,  certain  tidal  effects  can  also  be 
included  in  the  adjustment.  As  can  be  gathered  from[Blaha,  1982],  each 
diurnal  and  semidiurnal  tidal  constituent  considered  is  attributed  two 
adjustable  parameters,  a  global  amplitude  factor  and  a  global  phase  angle 
correction. 

The  altimeter  residuals  can  be  used  in  the  role  of  data  in  a  second- 
phase  solution  concerned  with  a  more  detailed  description  of  the  earth's 
gravity  field.  The  report  [Blaha,  1983]  makes  use  of  point-mass  parameters 
in  the  second-phase  adjustment  which  yields  a  banded  system  of  normal  equations 
resolved  through  the  application  of  the  "modified  Choleski  algorithm".  Due  to 
the  nature  of  this  algorithm,  the  geoid  over  entire  ocean  basins  can  be 
adjusted  in  a  few  overlapping  strips  of  point  masses,  leading  to  a  detailed 
resolution  of  the  oceanic  geoid  and  the  related  geophysical  quantities  on  a 
global  basis. 

The  material  contained  in  the  present  report  can  be  divided  into  four 
categories.  The  first  category  is  concerned  with  improvements  in  the  model 
of  the  sea  surface  and  of  tidal  effects;  the  second  category  is  concerned 
with  an  improvement  in  weights  attributed  to  altimeter  observations  in  the 
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first-phase  adjustment;  the  third  category  is  concerned  with  a  second-phase 
adjustment  which,  however,  includes  also  the  chosen  tidal  parameters  as 
adjustable  quantities;  and  the  fourth  category  is  concerned  with  an  applica¬ 
tion  of  the  least-squares  collocation  with  noise  to  a  second-  or  third-phase 
resolution  of  the  earth's  gravity  field. 

When  dealing  with  various  ocean-surface  characteristics  (first  category 
above),  one  recognizes  that  the  relationship  between  the  pressure  and  geo¬ 
potential  gradients  in  the  ocean  has  important  oceanographic,  meteorological 
and  geodetic  applications.  With  regard  to  the  usual  meteorological  practice, 
for  example,  such  a  relationship  leads  to  the  elimination  of  the  need  for  the 
variation  of  density  with  elevation  when  expressing  the  geostrophic  wind 
equations  on  an  isobaric  surface.  The  point  to  be  emphasized  here  is  that  the 
above  gradient  relationship  is  paramount  to  solving  important  geodetic  and 
oceanographic  problems,  such  as  the  determination  of  the  departures  of  the 
standard  ocean  from  the  geoid.  Without  taking  it  properly  into  account  one 
could  not  avoid  discrepancies  between  altimetric  and  oceanographic  solutions 
of  the  sea  surface.  The  refinement  of  the  usual  formula  giving  this  relation¬ 
ship  is  the  topic  of  Chapter  2. 

Since  the  altimeter  measurements  are  directly  affected  by  the  surface 
(geocentric)  tide,  an  exact  representation  of  the  latter  is  important.  The 
improvement  of  the  overall  tidal  model  in  its  relation  to  satellite  altimetry 
is  thus  linked  to  the  improvement  of  the  model  for  the  bottom  tide  which, 
when  added  to  the  ocean  tide,  yields  the  surface  tide.  This  topic  (also 
belonging  to  the  first  category)  is  addressed  in  Chapter  3. 

Due  to  an  imperfect  description  of  the  earth's  gravity  field  and  its 


fundamental  surface,  the  geoid,  by  a  truncated  spherical -harmonic  series 
in  a  first-phase  adjustment,  the  neglected  effect  is  superimposed  on  the 
observational  noise  and  acts  as  if  increasing  its  variance  and  introducing 
covariances  in  a  global  sense.  The  temptation  is  great  to  use  this  new 
variance,  and  to  neglect  the  geoidal  covariances.  The  advantage  achieved 
with  such  a  simplification  is  the  elimination  of  the  weight  matrix  of 
observations  from  the  least-squares  algorithm.  The  price  to  pay  for  this 
economical  gain  is  a  less  rigorous  adjustment.  A  compromise  solution  offer¬ 
ing  almost  equally  significant  computer  economies  is  presented  in  Chapter  4 
(complemented  by  Appendix  1). 

The  "modified  Choleski  algorithm"  developed  in  [Blaha,  1983]  is  based 
on  three  separate  steps:  1)  elimination  of  the  point  masses  from  an  observa¬ 
tion  equation  if  they  are  sufficiently  far  from  the  pertinent  observation 
point,  2)  special  arrangement  of  the  point-mass  parameters  in  the  adjustment 
scheme,  and  3)  resolution  of  the  resulting  system  through  an  adaptation  of 
the  well-known  Choleski  algorithm.  However,  this  approach  does  not  provide 
for  an  additional,  second-phase  adjustment  of  tidal  parameters  which  may  be 
desirable  at  the  level  of  individual  ocean  basins.  The  task  of  including 
chosen  tidal  parameters  in  the  second-phase  adjustment  of  point  masses  is 
addressed  in  Chapter  5. 

A  second-phase  adjustment  in  terms  of  parameters  such  as  the  point  masses 
requires  the  formation  of  normal  equations  and  their  resolution  in  a  simultane¬ 
ous  least-squares  process.  This  process  can  be  computationally  demanding, 
even  if  no  tidal  parameters  are  present  and  if  advantage  is  taken  of  the 
"modified  Choleski  algorithm".  An  equivalent  second-phase  approach  can  be 
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developed  in  using  an  adaptation  of  Moritz's  [1980]  least-squares  colloca¬ 
tion  with  noise.  This  approach  can  also  serve  in  predicting  other  geophy¬ 
sical  quantities  in  addition  to  geoid  undulations,  such  as  gravity  anomali 
It  is  described  in  Chapter  6  together  with  Appendices  2  and  3. 


2.  RELATIONSHIP  BETWEEN  THE  PRESSURE  AND  GEOPOTENTIAL  GRADIENTS 
IN  THE  OCEAN  DERIVED  IN  TENSOR  NOTATIONS 

The  increasing  accuracy  of  satellite  altimetry  and  its  growing  use 
in  the  determination  of  the  general  circulation  of  the  oceans  provide  the 
motivation  for  conceiving  the  most  rigorous  model  possible  in  relating 
the  pressure  and  geopotential  gradients.  It  is  deemed  worthwhile,  then, 
to  refine  the  standard  derivation  appearing  e.g.  on  pages  187-189  in 
[Hess,  1959]  by  including  what  could  be  termed  "second-order  effects". 

Even  if  these  effects  proved  insignificant  for  the  current  and  near-term 
prospects,  it  is  reassuring  to  know  the  extent  of  approximations  entering 
the  model  used  in  the  actual  applications. 

This  chapter  focuses  on  the  derivation  of  the  refined  gradient  formula 
in  curvilinear  coordinates  with  the  aid  of  tensor  notations.  Tt  then 
proceeds  to  the  simplification  removing  the  second-order  effects,  which 
results  in  the  standard  formula  as  presented  in  the  above  reference.  The 
differences  between  these  two  formulations  are  discussed  in  the  closing 
paragraphs. 

The  problem  to  be  addressed  is  depicted  in  Figure  1.  The  figure  displays 
the  curvilinear  nature  of  the  coordinates  associated  with  the  point  Q  (at 
the  intersection  of  the  two  planes  shown),  although  the  W-coordinate  line 
through  Q'  has  also  been  drawn.  In  solving  the  problem  at  hand,  some 
elements  of  tensor  analysis  as  presented  by  Hotine  [1969]  will  be  used. 

The  coordinate  system  employed  is  symbolized  by  xr,  r=l,2,3.  Here 
the  third  space  coordinate  is  x  =W,  the  geopotential;  x  and  x  are  surface 
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coordinates,  where  the  surface  is  W  =  constant,  as  well  as  the  other  two 
space  coordinates.  Accordingly,  xr=(x1,x2,W).  As  one  example,  x1  and  x2 
could  be  the  geodetic  latitude  and  longitude,  respectively. 

In  accordance  with  the  geodetic  convention,  the  potential  is  taken 
as  increasing  downward.  If  the  infinitesimal  vector  "dx"  associated  with 
the  point  Q  and  represented  by  its  contravariant  components  dxr  should  lie 
in  the  surface  W  =  constant,  then  along  its  direction  one  has 

dW  =  ( 9W/9xr)dxr  =  Wrdxr  =  0  . 

Thus  Wr  ,  the  gradient  of  W,  is  a  covariant  vector  perpendicular  to  the 
surface  W  =  constant.  If  represents  the  covariant  components  of  the 
(downward)  unit-vector  "v"  perpendicular  to  this  surface  at  Q,  the  above 
finding  is  formalized  as 

«r-9V  (2.1) 

where  g,  the  gravity  at  Q,  is  the  magnitude  of  . 

If  ds  is  the  length  along  v  so  that  the  contravariant  infinitesimal 
vector  considered  is  vrds  ,  then  along  this  direction  we  have 

dW  =  W  vrds  , 

r 

and  due  to  (2.1), 

dW  =  g  ds  ,  (2.2) 

which  is  the  familiar  geodynamic  equation. 
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W-coordlnate  line  (c.l.) 


Figure  1 

Schematic  representation  of  the  surfaces  VI  ■  constant  and  P  *  constant, 
where  W  Is  the  geopotentlal  and  P  Is  the  pressure  in  the  standard  ocean. 

The  coordinates  of  only  four  (out  of  six)  points  depicted  In  the  figure 
are  written  explicitly.  In  general,  the  zero  coordinates  could  be  re¬ 
placed  by  the  pertinent  x*  ,  x2  and  Wf0)  . 


So  far  the  coordinate  system  considered  is  quite  general  and  the 
direction  of  the  W-coordinate  line  (along  which  x  and  x  are  constant) 

does  not  have  to  coincide  with  v.  This  statement  applies  to  Q  or  any 

1  2 

other  point  in  the  surface  W  =  constant.  The  x  -  and  x  -  coordinate 
lines,  on  the  other  hand,  lie  in  this  surface.  Thus  the  W-coordinate 
line  is  not,  in  general,  perpendicular  to  x  -  or  x  -lines  at  Q  or  any 
other  point  of  the  surface. 

Consider  now  the  infinitesimal  neighborhood  of  Q  in  space.  Along 

1  2 

the  x  -,  x  -  and  W-coordinate  lines  we  have,  respectively: 

(2.3) 

dx*u  =  (dx1,0,0)  ,  dx^2)  =  (0,dx2,0)  ,  dx*3)  =  (0,0, dW). 

In  terms  of  Figure  1,  the  length  elements  of  these  infinitesimal  vectors 

2  r  s 

are  d£...  ,  i  =  1 ,2 ,3.  From  the  basic  relation  d£  =  g  dx  dx  giving  the 

( i )  rs 

square  of  the  length  element  along  any  vector  dx,  where  grg  is  the  metric 
tensor  at  Q,  one  finds 

^(1)  =  V  9ii"  dx  »  ^(2)  =v/922  ’  ^(3)  '^^33  ’ 

Since,  in  general,  the  contravariant  components  of  the  unit  vector  in  the 
direction  of  dx  are  £r  =  dxr/d£  ,  we  have  for  the  unit  vectors  i  in 
the  above  directions: 

(2.4) 

*u>  =  (WSIT.o.o)  •  *<2)  =  (O.IA^J.0)  .  »(3)  =  (0,0,1/^)  • 

These  formulas  are  valid  regardless  of  the  nature  of  the  coordinate  system. 
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Associated  with  the  interval  QQ*  of  Figure  1,  we  denote 

dxr  =  (dx1,dx2,0)  ,  (2.5) 

dZ  =  [g11(dx1)2  +  g22(dx2)2  +  2g12dx1dx2] 1/2  , 

=  (dx1/dA,dx2/d£,0)  .  (2.6) 

Since  xa,  a=l,2,  are  space  as  well  as  surface  coordinates,  it  holds  true 
that  g^  =  a^g  ,  where  aa^  is  the  surface  metric  tensor.  In  the  (x1,x2,W) 
coordinates  one  has  Wr  =  (0,0,1);  thus,  due  to  (2.1),  the  unit  vectors 
S.(1) ,  S,(2)  and  t  are  confirmed  to  lie  in  the  surface  W  =  constant. 

Next  introduce  normal  coordinates  so  that  the  pertinent  infinitesimal 
vectors,  which  are  here  split  naturally  according  to  their  coordinate 
differences  as 

(dx1,dx2,dW)  =  (dx\dx2,0)  +  (0,0, dW)  ,  (2.7) 

be  at  the  same  time  conveniently  split  into  a  vector  lying  in  the  surface 
W  =  constant  and  a  vector  perpendicular  to  it.  The  coordinates  x1  and  x2 
can  still  be  quite  general.  For  example,  they  could  again  be  the  geodetic 
latitude  and  longitude,  respectively,  but  they  would  have  such  characteristics 
only  on  one  surface,  in  this  case  the  above  surface  W  =  constant.  This 
transpires  from  the  paragraph  8  on  page  104  in  [Hotine,  1969],  where  such 
a  surface  is  called  "base  surface". 

1  2 

In  the  normal  coordinate  system  (x  ,x  ,W)  the  metric  tensor  can  be 
written  in  the  abbreviated  form  as 

9„  -  <9a8.l/92)  .  (2.8) 
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similar  to  equation  (15.03)  in  [Hotine,  1969].  This  means  that  in  the 
matrix  form  of  grs  the  elements  (1,3),  (2,3)  and  (3,1),  (3,2)  are  zero. 
Since  t  =  g  ts  holds  true  for  any  vector  t,  it  now  follows  from  (2.4) 

r  3rs 

that 

(2.9) 

*U>r’(*'9U’912^5I7*0)*  t<2lr'(912/'-'9I7>'^>0,>  *  ( 3)  r" 10'°' 1/9)  ’ 

where  g33  =  1/g  has  been  taken  into  account.  Upon  contractions  of  (2.9) 
with  (2.4)  and  (2.6),  it  is  verified  that  the  unit  vector  Jl  is 
perpendicular  to  the  unit  vectors  £(1)*  *-(2)  and  JL  We  thus  have 


vr  =  H(3)r  =  (0,0, l/g)  ,  vr  £  t*3)  =  (0,0, g)  , 


ds  =  ds,(3)  =  dW/g  . 


(2.10) 


(2.11) 


Upon  denoting  the  vector  along  QT'  as  dx' ,  where  T'  lies  in  the  surface 
P=constant,  this  vector  can  be  split,  according  to  (2.7),  as 


dx,r  =  dxr  +  dx^3)  , 


(2.12) 


where  dx'r  h  (dx1,dx2,dW)  corresponds  to  QT’  as  just  stated,  dxr  5  (dx1,dx2,0) 
corresponds  to  QQ'  (see  equation  2.5),  and  dx^3)  =  (0,0, dW)  corresponds  to 
QT  (see  equation  2.3).  According  to  the  foregoing,  dx  lies  in  the  surface 
W  --  constant  and  dx(3)  is  perpendicular  to  it.  All  three  vectors  dx',  dx 
and  dx(j),  whose  length  elements  in  Figure  1  are  dt' ,  d£  and  d5(3), 
respectively,  are  associated  with  Q. 

Consider  the  change  in  W  from  Q  to  T' ,  namely 


dW_,  =  W  dx,r  =  W  dxr  ♦  W  dxr 

QT'  r  r  r  (  i) 


-10- 


The  first  term  on  the  right-hand  side  is  zero  due  to  the  relation  preceding 
(2.1).  In  consulting  (2.2)  and  the  equation  that  preceded  it,  the  second 
term  with  dx*3)  =  vrds  =  t^3)d«,(3)  leads  to 
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(2.13) 


Next  consider  the  change  in  P  from  Q  to  Q'.  Due  to 


dP  =  P  dx,r  =  P  dxr  +  P  dx*  =  0  , 

QT*  r  r  r  (3) 

it  follows  that 


P  dxr  =  -P  dx* 
r  r  (3) 


However,  the  hydrostatic  equation  for  the  compressible  standard  ocean  at  Q 
yields 


-  9  0  "(3,  •  <2'I4> 

where  p  is  the  water  density  at  Q,  due  to  the  fact  that  dJl(3)  is  the  dis¬ 
placement  in  the  direction  of  gravity  (i.e.,  perpendicular  to  the  equi- 
potential  surface  W  =  constant).  Accordingly, 

dV  ’  -9  C  dll<3,  •  <2-I5> 

In  combining  (2.13)  and  (2.15),  we  have 


dWQT,  =  -(l/pJdPgg.  ,  (2.16) 

specifying  the  change  in  geopotential  along  an  isobaric  surface  in  terms 
of  the  change  in  pressure  along  an  equipotential  surface.  Next  we  address 
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the  problem  of  the  rate  of  change  in  these  quantities  along  the 
respective  paths. 

The  following  relations  are  obtained  from  (2.13)  and  (2.15), 
respectively: 


M/di'  =  g(dt(3)/ds.' )  , 


3P/3*.  =  -pg(d«,(3)/dJt)  . 


(2.17) 


(2.18) 


These  equations  yield 


3W/3H'  =  -(1/p)  f  3P/3£  , 


f  =  dil/dSL'  . 


(2.19) 


From  the  formula  d£'  =  9rsdx'  dx's  ,  where  g^g  is  given  in  (2.8) 

and  dx‘r  is  given  following  (2.12),  one  obtains 


dll'  =  CdH2  +  (dW/g) 2  ]1/2  s  (dll2  +  d?23))1/2  , 


(2.20) 


where  d£  appears  in  the  relations  following  (2.5)  and  dH(3)  appears  in  (2.11), 


From  (2.19)  and  (2.20),  one  writes 


f  =  (1  +  dC(3)/dty 


[1  +  (dW/g)2/dsi2;r 


(2.21) 


If  the  two  surfaces  are  "close"  to  each  other,  in  the  sense  that 


dn(3)  h  dW/g  «  dH  , 


then  the  approximations 


f  =  1  , 

aw/ar  =  -(l/pjap/at 


(2.22) 


✓ 

are  satisfactory.  In  this  context  one  can  regard  the  pressure  gradient 
along  equi potential  surfaces  as  interchangeable  with  the  geopotential 
gradient  along  isobaric  surfaces  (except  for  a  factor),  a  classical 
oceanographic  procedure. 

When  comparing  equation  (2.22)  with  the  corresponding  equations  in 
[Hess,  1959]  (it  precedes  equation  12.10  therein  and  will  be  designated 
by  the  symbol  "*"),  we  notice  the  difference  in  sign.  However,  this 
stems  from  the  orientation  of  the  z-axis  in  *  (positive  upward)  as  opposed 
to  the  distance  d£(3)  taken  here  as  positive  downward.  If,  in  the  above 
reference,  one  stipulated  that  dW  =  -gdz  for  W  increasing  downward  in 
agreement  with  the  geodetic  convention,  the  difference  in  sign  would 
disappear. 

The  equivalent  of  *  can  be  obtained  upon  dividing  both  sides  of 
equation  (2.16)  by  d£.  The  left-hand  side  of  this  new  equation  becomes 
the  geopotential  gradient  along  an  isobaric  surface,  9W/9J,'  ,  if  we  assume 
that  dn=d£'.  With  this  approximation  the  new  equation  is  then  (2.22). 
Clearly,  a  simplification  of  this  kind  makes  the  rigorous  formulation 
(2.19)  and  the  practical  formulation  (2.22)  indistinguishable.  Thus  the 
standard  derivation  leading  to  *  causes  the  rigorous  gradient  model  (2.19) 
with  f  given  by  (2.21)  to  remain  hidden.  As  we  have  seen,  such  a  relation¬ 
ship  can  be  unearthed  if  one  considers  the  gravity- related  (curvilinear) 
coordinate  system  and  focuses  on  its  W-coordinate. 
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3.  IMPROVED  FORMULATION  OF  THE  SURFACE  TIDE  IN  THE  SATELLITE 

ALTIMETRY  ADJUSTMENT  MODEL 

The  altimeter  measurements  are  directly  affected  by  the  sufface 
(geocentric)  tide.  Since  the  surface  tide  is  given  by  the  addition  of  the 
ocean  tide  and  the  bottom  tide,  an  exact  representation  of  the  latter  is 
important.  However,  when  the  tidal  adjustment  was  formulated  in  [Blaha,  1982], 
certain  approximations  were  adopted,  such  as 

-u^  =  0.06£  .  M2  constituent  , 

-u  =0  .  all  the  other  diurnal  and 

semidiurnal  constituents  , 

where  u^,  a  part  of  the  bottom  tide,  is  the  ocean  bottom  deformation  due  to  the 
ocean  tidal  loading  and  £  is  the  ocean  tide.  These  approximations  were  made 
through  the  inspections  and  comparisons  of  tidal  contour  maps,  especially 
those  pertaining  to  the  M2  tide  (see  pages  97-106  of  [Estes,  1980]  and  pages 
386-402  of  [Parke  and  Hendershott,  1980]). 

In  using  (3.1),  the  following  model  was  formulated  on  page  56  of  [Blaha, 
1982]  for  the  diurnal  and  semidiurnal  constituents  identified  by  the  index  j: 

Ik  =  c'c^  +  h  x  (equilibrium  tide)..  ,  (3.2a) 

where 

c'  =  0.94  .  M2  constituent, 

c'  =  1.00  .  the  other  diurnal  and 

semidiurnal  constituents, 

h  =  0.61,  a  Love  number. 
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and  where 


h..  =  surface  tide  for  the  constituent  j  , 

5.  =  ocean  tide  for  the  constituent  j  . 

D 

The  surface  tide  was  considered  in  the  above  reference  on  pages  39-42,  where 
it  was  denoted  as  5  . 

s 

An  improvement  in  the  above  model  can  be  achieved  by  adopting  a  more 
complex  tidal  model  developed  by  Schwiderski  [1980],  where 

5b  =  5e  +  5*°  ,  (3.3) 

wi  th 

£b  =  bottom  tide, 

5e  =  earth  tidal  elevation  expressed  as  0.61n,  n  being  the  Newton's 
equilibrium  tide  and  0.61  being  the  Love  number  h, 

5eo=  earth-dip  response  to  oceanic  tidal  load. 

The  tidal  model  has  been  further  refined  in  Schwiderski *s  recent  NSWC  reports 
concerned  with  the  global  ocean  tides.  Although  these  reports  do  not  deal 
with  the  altimeter  adjustment  per  se,  certain  features  of  the  tidal  model 
presented  therein  can  have  a  beneficial  effect  on  the  global  adjustment  model 
of  satellite  altimetry. 

In  particular,  the  following  improved  relationships  developed  by 
Schwiderski  can  be  used  with  advantage: 

5*  =  0.612n  ,  (3.4a) 

5eo  =  -0.06675.  (3.4b) 
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Equation  (3.4b)  provides  an  improved  formula  for  the  ocean  bottom  deformation 
due  to  ocean  tidal  loading,  written  as 

-ul  =  0.0667S  ,  (3.4b') 

which  can  be  used  in  conjunction  with  all  of  the  tidal  constituents.  Combined 
with  (3.3),  equations  (3.4a,b)  imply  that 

£b  =  0.612n  -  0.0667C  .  (3.5) 

Due  to 

Cs  =  £  +  Kh  , 

where  Cs  is  the  surface  tide  (in  [Blaha,  1982]  s  and  b  were  used  as  lower 
indices),  it  follows  that 

Cs  =  0.9333C  +  0.612n  .  (3.6) 

The  final  formula  (3.6)  can  be  used  in  conjunction  with  all  the  constituents, 
long-period  as  well  as  diurnal  and  semidiurnal.  In  the  notations  of  [Blaha, 

1982],  it  can  be  transcribed  as 

hj  =  c'Cj  +  h  x  (equilibrium  tide)^  ,  (3.7a) 

which  has  the  same  form  as  (3.2a),  but  where 

c'  =  0.9333  .  all  of  the  tidal  constituents.  j 

>  (3.7b) 

h  =  0.612,  a  Love  number. 
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The  improvement  achieved  with  the  refined  model  can  be  assessed  by  comparing 
(3.7b)  with  (3.2b).  The  change  in  the  value  of  h  is  of  little  consequence. 
However,  the  change  affecting  the  M2  constituent  is  more  important  as  are 
the  changes  affecting  the  other  diurnal  and  semidiurnal  constituents  (the 
relative  changes  affecting  the  latter  are  greater  than  the  relative  change  in 
M2  ,  but  this  relationship  does  not  necessarily  hold  true  for  the  absolute 
changes  in  the  geocentric  distance).  Due  to  the  changes  in  c'  above,  the 
partial  derivatives  with  respect  to  all  the  diurnal  and  semidiurnal  constitu¬ 
ents  are  modified  as  well. 

The  long-period  constituents  have  been  treated  in  [Blaha,  1981,  1982] 
in  a  simplified  model,  where  the  ocean  tidal  loading  as  well  as  other  effects 
were  disregarded.  The  model  was  based  on  a  modified  equilibrium  tide  since 
no  spherical -harmonic  tidal  coefficients  were  available  to  express  the  ocean 
tide  £;  for  these  constituents.  However,  this  situation  may  change  due  to 
Schwiderski 's  recent  reports  offering  a  possibility  to  derive  such  coef¬ 
ficients  and  thus  to  treat  some  or  all  of  the  long-period  constituents  Mf, 

Mm  and  SSa  in  the  same  fashion  as  the  diurnal  and  semidiurnal  constituents 
above,  i.e.,  via  equations  (3.7a,b).  We  can  close  this  discussion  by  stating 
that  the  modifications  (3.7b)  represent  an  improvement  of  the  tidal  adjustment 
model  expected  to  be  translated  into  an  improvement  in  satellite  altimetry 
results  by  several  centimeters,  especially  in  the  tidal  effects  themselves 
and  in  the  altimetry  residuals. 
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4.  WEIGHTING  SCHEME  FOR  A  GLOBAL  ADJUSTMENT 
OF  SATELLITE  ALTIMETRY 


The  data  in  satellite  altimetry  consist  of  altimeter  observations  whose 
noise  level  has  been  reduced,  in  the  case  of  SEASAT,  to  0.1-0. 2m.  The  global 
(first-phase)  adjustment  model  provides  for  a  simultaneous  least-squares 
solution  of  a  truncated  set  of  spherical -harmonic  (S.H.)  potential  coef¬ 
ficients,  the  state  vector  (s.v.)  parameters  on  all  satellite  arcs  and, 
optionally,  certain  tidal  parameters.  However,  due  to  an  imperfect  description 
of  the  earth's  gravity  field  and  its  fundamental  surface,  the  geoid,  by  a 
truncated  series,  the  neglected  effect  will  be  superimposed  on  the  observa¬ 
tional  noise  and  act  as  if  increasing  its  variance  and  introducing  covariances 
in  a  global  sense.  Therefore,  if  the  altimeter  observations  were  treated  as 
independent  and  as  having  the  one-sigma  level  of  0.1-0. 2m,  they  would  largely 
override  any  conflicting  information  supplied  by  the  weighted  S.H.  coefficients 
and,  especially,  the  weighted  s.v.  parameters.  The  problem  would  be  further 
compounded  by  high  density  of  geoidal  data. 

The  last  statement  can  be  corroborated  by  the  realization  that  the  geoid 
represents  a  smooth  function  and  the  participation  of  a  great  number  of 
observations  short  distances  apart  would  lead  to  a  similar  overall  effect  on 
the  above  parameters  (the  observations  treated  as  independent  would  be  allowed 
to  unduly  reinforce  each  other).  This  effect  can  be  inferred  from  the  geoidal 
covariance  function  constructed  for  a  given  truncation  level,  allowing 
quantitative  evaluations  of  both  the  geoidal  variances  at  single  points  and 
the  geoidal  covariances  for  pairs  of  points  separated  by  given  distances. 

The  lower  the  degree  and  order  truncation,  the  larger  these  numbers  and  the 
smoother  the  covariance  function. 
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The  first-phase  adjustment  at  AFGL  is  performed  in  terms  of  a  (14,14) 
truncated  set  of  S.H.  potential  coefficients.  The  covariance  function  can 
be  used  here  to  quantify  the  overall  effect  of  the  neglected  coefficients, 
from  degree  n=  15  onward  (in  theory,  to  infinity;  in  practice,  to  n=  500 
or  n=  1,000,  etc.).  This  effect  can  be  represented  by  the  one-sigma  level, 
i.e.,  the  square-root  of  the  geoidal  variance,  which  for  the  (14,14)  trunca¬ 
tion  has  been  computed  as  4.9m.  However,  in  using  SEASAT  observations  this 
quantity  is  estimated  as  3.2m;  the  procedure  giving  this  estimate  will  be 
presented  later.  Combining  this  value  with  the  0.1-0. 2m  sigma  of  the  observa¬ 
tional  noise  proper  still  yields  about  3.2m.  Even  though  the  total  input 
sigma  is  thus  made  substantially  higher  than  0.1-0. 2m,  the  neglect  of  the 
geoidal  covariances  could  still  be  detrimental  to  the  adjustment  in  the  above 
sense  of  undue  reinforcement  of  observations.  But  the  temptation  to  neglect 
the  geoidal  covariances  is  great  from  the  practical  standpoint.  The 
significant  advantage  achieved  with  such  a  simplfi cation  is  the  elimination 
of  the  weight  matrix  of  observations  from  the  least-squares  algorithm  (more 
precisely,  the  reduction  of  a  fully  populated  weight  matrix  for  each  satellite 
arc  to  a  constant  times  the  unit  matrix).  A  compromise  solution  offering 
almost  equally  significant  computer  economies  will  be  presented  in  this 
chapter. 

The  proper  treatment  of  altimeter  observations  would  consist  in  filling 
the  input  variance-covariance  matrix  with  the  values  obtained  through  the  use 
of  the  geoidal  covariance  function  and  termed  "rigorous".  In  choosing  this 
avenue,  one  would  be  faced  with  the  task  of  inverting  one  such  matrix  per 
satellite  arc,  or  pre-inverting  such  matrices  according  to  the  number  of 
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observations.  If  all  the  arcs  had  the  same  number  of  observations  —  such 
as  54  at  0.5°  intervals  making  up  for  about  26°  arcs  at  seven-minute  durations 
which  has  been  the  upper  limit  adopted  for  SEASAT  —  only  one  pre-inverted 
matrix  (the  weight  matrix  P)  would  be  needed.  However,  this  is  not  the  case 
in  practice.  The  subsequent  formation  of  normal  equations  with  the  fully 
populated  weight  matrix  would  represent  further  requirements  compared  with  the 
simplified  approach  of  neglected  covariances.  For  example,  a  greatly  increased 
computer  storage  would  be  needed  because  the  whole  design  matrix  (A)  of 
observation  equations  for  an  arc  would  have  to  be  retained  in  the  active 
memory.  Its  dimensions  are  (m*u),  m  being  the  number  of  observations  and  u 
being  the  number  of  parameters.  In  the  simplified  approach  only  one  row  of 
this  matrix  is  needed  at  a  time.  With  regard  to  the  computer  run-time  during 
the  formation  of  the  matrix  of  normal  equations  associated  with  that  arc, 
the  number  of  scalar  multiplications  would  also  increase  significantly  in 
such  a  comparison,  from  mu2  in  the  simplified  approach  to  mu2+m2u  in  the 
"rigorous"  approach. 
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4.1  First-Order  Autoregressive  Process 


A  compromise  solution  will  now  be  addressed,  promising  from  both  the 
accuracy  and  the  economy  standpoints.  The  input  variance-covariances  will  be 
slightly  modified  from  their  "rigorous"  values  as  represented  by  the  geoidal 
covariance  function.  As  a  result,  a  highly  patterned  tri -diagonal  weight 
matrix  known  beforehand  will  become  available.  This  "modified  approach" 
corresponds  essentially  to  the  first-order  autoregressive  process  described 
in  [Brown  and  Trotter,  1969].  In  the  same  reference,  the  second-order  auto¬ 
regressive  process  would  correspond  to  a  five-diagonal  weight  matrix,  etc. 

In  the  case  at  hand,  the  tri-diagonal  weight  matrix  is  arrived  at  via  slight 
deformations  of  the  original  covariance  function.  In  particular,  if  the 
covariance  function  becomes  of  the  form  yielding  the  following  variance- 
covariance  matrix. 


C  =  a 


1 


„  r?  _3  m-1 

P  P  P  .  P 

„  i  „  „2  m-  2 

P  1  P  P  .  P 

2  ,  m-3 

P  P  1  P  .  P 

m-1  m-2  , 

P  P  .  1 


(4.1) 
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then  the  weight  matrix  follows  as 


1  -p  0 

-P  1+p2  -p 


P  H  C'1  =  (l/[o2(l-p2)]} 


0  -p  1+p2 


0 

0 

-p 


(4.2) 


0  0  .  -P  1+p2  -p 

0  0  .  0  -p  l 


as  can  be  verified  upon  pre-  or  post-multiplying  C  by  P. 

Clearly,  when  passing  from  the  diagonal  (simplified)  to  the  tri-diagonal 
(modified)  case,  the  increases  in  the  storage  requirements  and  in  the  number 
of  scalar  multiplications  are  negligible  when  compared  to  the  increases 
associated  with  the  "rigorous"  case.  Furthermore,  the  pattern  in  (4.2)  is 
the  same  regardless  of  the  dimensions  of  P  so  that  no  separate  inversions,  or 
storage  of  pre-inverted  matrices,  are  necessary  for  m  varying  from  arc  to  arc, 
unlike  in  the  "rigorous"  case.  Yet,  the  tri-diagonal  case  yields  the  solution 
which  is  likely  to  closely  approximate  the  "rigorous"  solution,  much  more  so 
than  could  be  expected  from  the  diagonal  case.  This  depends,  of  course,  on 
how  much  the  covariance  function  is  to  be  deformed  in  order  to  fit  the  tri¬ 
diagonal  scheme.  A  simple  procedure  of  arriving  at  the  tri-diagonal  case  is 
presented  in  the  sections  below.  The  analysis  will  be  concluded  with  a 
special  example  showing  the  improvement  represented  by  the  tri-diagonal  case 
versus  the  diagonal  case. 
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4.2  Geoidal  Covariance  Function 


The  representation  of  the  covariance  function  for  geoid  undulations  is 
based  on  the  results  of  [Tscherning  and  Rapp,  1974].  These  results  were 
adapted  in  [Blaha,  1982]  to  give  the  k-degree  variance  (in  meters  square)  as 
fol 1 ows : 

o2  =  17,981  x  0.999617k+2/ [(k-l)(k-2)(k+24)] .  (4.3) 

The  covariance  function  for  geoid  undulations  is  then  given  as 

00 

D(<|>)  =  M{N,N'  }  =  l  a2  P.  (cosi/O  , 
k=3  *  * 

where  the  operator  M  indicates  averaging  over  the  whole  unit  sphere,  N  and  N’ 
are  the  geoid  undulations  at  any  two  points  P  and  P'  separated  by  the  spherical 
distance  i|>,  and  Pk(cosij;)  are  the  Legendre  polynomials  in  the  argument  cos4<. 

The  covariance  function  corresponding  to  the  bandwidth  n1  through  n2  would 
read 

02  2 

O'M  =  I  \  Pk(cosif;)  , 

k=n^ 

whereas  the  covariance  function  corresponding  to  the  truncation  n  can  be 
symbolized  by 

00 

cov(n+l,ij;)  =  l  a2  Pk(cos^)  .  (4.4a) 

k=n+l 

Here  "cov"  depicts  the  geoidal  covariances,  while  "var"  depicts  the  geoidal 
variances,  both  due  to  the  neglected  degrees  (from  n+1  to  infinity).  Since 
Pk(l)  =  1.  one  has 
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(4.4b) 


var(n+l)  =  l  a2  . 

k=n+l 

In  practical  computations  the  degree  is  substituted  for  by  n^^  which, 
as  indicated  earlier,  can  be  replaced  by  500  or  1,000,  or  some  other 
sufficiently  large  number. 

Appendix  1  features  a  computer  program  using  the  algorithm  (4.4a,b), 
where  the  degree  variances  are  given  in  (4.3)  and  the  well-known  recurrent 
relationship  for  the  Legendre  polynomials  reads 

Pk(z)  =  (l/k)[z(2k-l)Pk_1(z)  -  (k-l)Pk_2(z)]  , 

starting  with 

P0(z)  =  1  ,  Px(z)  =  z  . 

This  formula  is  contained  in  the  subroutine  "Legen"  also  listed. 

The  symbols  utilized  by  the  program  in  Appendix  1  are: 

NPROB  ...  number  of  different  problems  to  be  treated, 
and  within  each  problem: 

IORDER  ...  n  (dimensioned  for  1,000), 

max 

NPSI  ...  number  of  different  angles  V  (dimensioned  for  30), 

NN  ...  number  of  different  truncations  "n"  (dimensioned  for  20). 

Following  NPSI  the  input  consists  of  the  individual  angles  In  degrees,  and 
following  NN  the  input  consists  of  the  individual  truncations  n  in  the 
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descending  order.  In  the  computer  run  carried  out  for  the  present  purpose 
NPROB  has  been  1,  I ORDER  has  been  500,  NPSI  has  been  27,  the  (27)  values  of 
ij;  have  run  from  0  to  13°  in  0.5°  intervals,  NN  has  been  1,  and  the  (one) 
desired  truncation  has  been  n=  14.  For  comparison  purposes,  the  values 
associated  with  ij>  =  15°  and  \p=  30°  will  also  be  presented. 

The  values  corresponding  to  the  (14,14)  truncated  model  as  obtained  with 


the  above  computer  program  are  (t{;  is  in 

°,  COV\p 

is  in  nr 

«#  • 

-h 

O 

-5 

II 

0  "cov" 

means 

"var" ) : 

(4.5) 

: 

0 

0.5 

1.0 

1.5 

2.0 

2.5 

3.0 

3.5 

4.0 

cov\p: 

23.869 

23.102 

21.496 

19.501 

17.315 

15.055 

12.798 

10.600 

8.503 

<P  : 

4.5 

5.0 

5.5 

6.0 

6.5 

7.0 

7.5 

8.0 

8.5 

covijj: 

6.533 

4.710 

3.050 

1.560 

0.244 

-0.897 

-1.866 

-2.667 

-3.306 

: 

9.0 

9.5 

10.0 

10.5 

11.0 

11.5 

12.0 

12.5 

13.0 

covif>: 

-3.792 

-4.136 

-4.347 

-4.437 

-4.418 

-4.304 

-4.105 

-3.835 

-3.507 

The  covariances 

whose  magnitude 

is  less 

than  1/3  of  the 

variance 

(23.869m2) 

are  deemed  relatively  unimportant  for  the  adjustment.  Thus  the  factor  p 
appearing  in  the  modified  matrix  C  in  (4.1)  will  be  computed  using  only  the 
first  nine  values  (for  ranging  from  0  to  4°).  The  covariance  values  at 
^=4°  and  ^=4.5°  correspond  to  35%  and  27%  of  the  variance,  respectively.  The 
covariance  value  at  iJf130  does  not  even  attain  15%  of  the  variance,  and  the 
covariance  magnitudes  diminish  quite  rapidly  beyond  this  angle.  For  example, 
for  ^15°  the  covariance  value  Is  -1.841  (under  8%  of  the  variance)  and  for 
^=30°  this  value  Is  -0.574  (about  2%  of  the  variance).  The  values  of  the 
modified  covariance  function  corresponding  to  (4.1)  would  approach  zero  even 
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faster.  It  thus  appears  to  be  amply  sufficient  to  extend  this  analysis  as 
far  as  ^=13°  and  not  beyond. 


4.3  Estimated  and  Modified  Variance-Covariances 


As  has  been  already  indicated  in  [Blaha,  1982],  the  geoidal  variance 
from  (4.5)  --  or,  equivalently,  its  square  root  o=  4.886m  —  appears  to  be 
conservative  when  compared  to  its  estimate  using  SEASAT  altimeter  observa¬ 
tions.  The  final  root  mean  square  (rms)  of  the  constant  terms  of  observa¬ 
tion  equations  from  all  the  SEASAT  arcs  in  a  (14,14)  S.H.  model  has  been 
computed  as  3.6064m  and  designated  as  the  estimate  o(total).  The  estimate 
of  the  variance  (a2)  due  to  the  truncation  is  then  found  as 

o2  =  o2(total)  -  o2('.v.)  -  o2(sea  surf.)  -  o2(altim. )  -  o2(S.H.)  -  a2(algor.), 

where  o(s.v.)  is  the  estimated  sigma  due  to  the  s.v.  parameters,  in  particular, 
due  to  the  uncertainty  in  the  "up"  component  of  the  precise  ephemeris  used 
and  given  a  priori  as  1.6m;  a(sea  surf.)  is  the  estimated  sigma  due  to  the 
tides  and  other  sea  surface  effects,  assumed  to  be  in  the  vicinity  of  0.3m; 
o(altim.)  is  the  estimated  sigma  of  the  altimeter  noise,  assumed  to  be  0.2m; 
o(S.H.)  is  the  estimated  sigma  due  to  the  errors  in  the  a  priori  values  of  the 
(14,14)  set  of  S.H.  potential  coefficients  used,  whose  value  is  unknown  but 
considered  small  and  adopted  as  zero;  and  a(algor. )  is  the  estimated  sigma 
due  to  the  short-arc  algorithm,  also  considered  small  and  adopted  as  zero. 

With  the  above  values  one  finds 

a2  =  10.3161m2  ,  (4.6) 

which  is  substantially  smaller  than 

a2  =  23.869m2  ,  (4.7) 
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as  given  in  (4.5).  From  (4.6)  one  has  a=  3.212m,  which  may  still  be  slightly 
conservative  due  especially  to  the  assumption  o(S.H.)=  0.  In  any  event,  the 
effect  of  o(s.v.)  represents  the  largest  reduction  of  a(total);  with  only 
these  two  quantities  present  the  value  of  a  would  have  been  3.232m,  differing 
little  from  3.212m  computed  above. 

Due  to  the  lack  of  other  indications,  we  assume  that  the  variance- 

*2  2 

covariances  of  (4.5)  should  be  scaled  down  by  the  factor  o  /o  computed  from 
(4.6)  and  (4.7)  in  order  to  yield  what  will  be  called  the  estimated  geoidal 
variance-covariances.  This  factor  is 

a2 /a2  =  0.4322.  (4.8) 

/v 

The  ratio  a/a  is  0.657,  or  the  estimated  sigma  is  about  2/3  of  the  theoretical 
sigma.  In  applying  the  factor  in  (4.8)  to  (4.5)  one  obtains  (4.9)  below  where, 
however,  the  estimated  variance  of  10.3161m  has  been  increased  by  (0.4m)  . 

This  increase  is  due  to  the  altimeter  noise  represented  by  the  variance 
(0.2m)  ,  as  well  as  to  other  modeled  and  unmodeled  effects.  The  sigma  cor¬ 
responding  to  this  final  variance  is  3.237m.  The  resulting  estimates  are 

(4.9) 


<1*  : 

0 

0.5 

1.0 

1.5 

2.0 

2.5 

3.0 

3.5 

4.0 

covtp: 

10.476 

9.985 

9.291 

8.428 

7.484 

6.507 

5.531 

4.581 

3.675 

ip  : 

4.5 

5.0 

5.5 

6.0 

6.5 

7.0 

7.5 

8.0 

8.5 

covtp: 

2.824 

2.036 

1.318 

0.674 

0.105 

-0.388 

-0.806 

-1.153 

-1.429 

tp  : 

9.0 

9.5 

10.0 

10.5 

11.0 

11.5 

12.0 

12.5 

13.0 

covty: 

-1.639 

-1.788 

-1.879 

-1.918 

-1.909 

-1.860 

-1.774 

-1.657 

-1.516 
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As  indicated  earlier,  only  the  values  between  i^=0  and  <j/=4°  will  be  used 
in  computing  the  modified  variance-covariances  from  the  estimated  variance- 
covariances  of  (4.9).  In  particular,  a  coefficient  p  is  sought  which  gives 
a  covariance  in  (4.9)  upon  the  multiplication  by  the  preceding  value.  The 
eight  successive  ratios  obtained  from  the  first  nine  consecutive  values  in 
(4.9)  are 

0.9531,  0.9305,  0.9071,  0.8880,  0.8695,  0.8500,  0.8282,  0.8022. 

Beyond  4°  the  ratios  would  decrease  more  rapidly  (0.7684,  0.7210,  0.6473,  etc). 
The  average  coefficient  computed  from  the  above  eight  values,  and  its  second 
power,  are 

p  =  0.8786  ,  p2  =  0.7719  .  (4.10a,b) 

Upon  applying  (4.10a)  to  10.476  in  (4.9),  the  modified  covariances  are,  in 
succession, 

(4.11) 


^  : 

0 

0.5 

1.0 

1.5 

2.0 

2.5 

3.0 

3.5 

4.0 

covtp: 

10.476 

9.204 

8.087 

7.105 

6.243 

5.485 

4.819 

4.234 

3.720 

it  : 

4.5 

5.0 

5.5 

6.0 

6.5 

7.0 

7.5 

8.0 

8.5 

covtj >: 

3.268 

2.871 

2.523 

2.217 

1.948 

1.711 

1.503 

1.321 

1.161 

i>  : 

9.0 

9.5 

10.0 

10.5 

11.0 

11.5 

12.0 

12.5 

13.0 

co  vir. 

1.020 

0.896 

0.787 

0.692 

0.608 

0.534 

0.469 

0.412 

0.362 
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At  this  point  three  weighting  schemes  for  observational  intervals  of  0.5°, 
sufficient  for  a  1°-  or  a  2°-geoidal  resolution,  are  presented.  The  first, 
"rigorous"  case  is  characterized  by  the  C  matrix  of  dimensions  (mxm),  where 
m  is  now  taken  as  27,  whose  entries  correspond  to  the  estimated  variance- 
covariances  of  (4.9): 


— 

10.476 

9.985 

9.291  ... 

-1.657 

1 - 

vo 

H 

LT> 

T—i 

1 

9.985 

10.476 

9.985  ... 

-1.774 

-1.657 

9.291 

• 

9.985 

10.476  ... 

-1.860 

-1.774 

• 

• 

-1.516 

-1.657 

-1.774  ... 

9.985 

10.476 

— 

The  inversion  of  this  matrix  yields  a  fully  populated  P  matrix. 

The  second,  tri -diagonal  case  is  characterized  by  the  C  matrix  of 
dimensions  (mxm),  where  again  m=  27,  whose  entries  correspond  to  the  modified 
variance-covariances  of  (4.11),  namely 


10.476 

9.204 

8. 087  ... 

0.412 

0.362 

9.204 

10.476 

9.204  ... 

0.469 

0.412 

8.087 

9.204 

10.476  ... 

0.534 

0.469 

0.362 

0.412 

0.469  ... 

9.204 

10.476 
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The  tri-diagonal  weight  matrix  is  given  as  in  (4.2),  where  the  value  "o2" 
is  10.476  and  p,  p2  are  given  in  (4.10a,b).  The  result  is 


P2  = 


a  b 
b  c 
0  b 


0  0  0  ...  0  0 


b  0  0 
c  b  0 


,..  0  0 

,..  0  0 


0  0  0  .  b  a 


(4.12) 


with 


a  =  0.4185,  b  =  -0.3677,  c  =  0.7415 


(4.13) 


The  third,  diagonal  case  is  characterized  by  the  simplification  in 
which  all  the  covariances  are  ignored,  resulting  in  the  following  matrices: 

C3  =  10.476  I  , 

P3  =  0.09546  I  ,  (4.14) 

where  I  is  the  identity  matrix  of  dimensions  (mxm),  again  with  m=  27. 

In  the  first  case,  the  formation  of  normal  equations  for  one  arc  con¬ 
taining  m  observations  (on  SEASAT  arcs  m?54)  proceeds  without  any  simplifi¬ 
cation  according  to  the  standard  formula 

(A^Ajx  =  A^e  ,  (4.15) 
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m 


where  A,  consisting  of  the  rows  A^  i=  1,2,. ..m  ,  is  the  (design)  matrix  of 
observation  equations  and  e  is  the  vector  containing  m  constant  terms  (ei), 
x  being  the  vector  of  parameters.  Here  all  the  m  rows  of  A  have  to  be  re¬ 
tained  in  the  computer  core  as  has  already  been  indicated.  The  same  holds 
true  for  the  m  constant  terms. 

In  the  second  case,  the  matrix  of  normal  equations,  ATP2A,  is  formulated 
in  a  much  simpler  manner,  namely 

ATP  A  =  (aA^  bA^)A  +mj1  [cAT+b(AT  +  AT  )]A.  +  (aAT+  bAT  )A  .  (4.16) 

2  1  2  1  .  *•_  l  l-l  i+l  l  m  m-1  m 

1=2 

T 

The  right-hand  side  of  normal  equations,  A  P2e,  is  computed  as  above,  except 
that  each  symbol  A  ,  k=  l,2,...m,  is  replaced  by  the  corresponding  e  .  The 

K  K 

middle  terms  on  the  right-hand  side  of  (4.16)  indicate  that  only  three  rows 
of  A  and  the  corresponding  three  elements  of  e  are  needed  in  the  computer  core 
at  any  one  time.  Clearly,  the  number  of  scalar  multiplications  in  this  process 
is  greatly  reduced  when  compared  to  the  first  case. 

The  third  case  is  the  simplest,  where 

ATP  A  =  0.09546  ?  ATA.  .  (4.17) 

3  .‘•.11 

1=1 

The  computation  of  ATP3c  proceeds  as  above  except  that  each  symbol  is  to 
be  replaced  by  the  corresponding  ei-  Here  only  one  row  of  A  and  the  cor¬ 
responding  element  of  e  are  needed  in  the  computer  core  at  any  one  time.  The 
number  of  scalar  multiplications  is  lower  than  in  the  previous  case,  but  the 
difference  is  not  great.  In  fact,  the  decrease  is  approximately  equal  to  the 
(2/u)-multiple  of  the  total  number  of  scalar  multiplications,  where  u  is  the 
number  of  parameters. 
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4  —.'Hm  wmi  *  •  ■- 


Although  the  second  case  is  about  as  economical  as  the  third  case, 
its  results  are  likely  to  be  much  closer  to  the  "rigorous"  (first)  case  than 
the  results  of  the  third  case.  At  first  sight,  this  can  be  assessed  as 


follows.  Let  the  differences  between  the  second  and  first  cases  be  called 
the  second  case  errors.  Upon  comparing  the  values  in  (4.11)  with  those  in 
(4.9),  these  errors  (between  ^=0  and  ^=13°)  are 


0, 

-0.781, 

-1.204, 

-1.323, 

-1.241, 

-1.022, 

-0.712, 

-0.347, 

(4.18) 

0.045, 

0.444, 

0.835, 

1.205, 

1.543, 

1.843, 

2.099, 

2.309, 

2.474, 

2.590, 

2.659, 

2.684, 

2.666, 

2.610, 

2.517, 

2.394, 

2.243, 

2.069, 

1.878. 

For  15° 

and  30° 

these  values  would 

be  1.012 

and  0.252 

,  respectively. 

The  rms 

of  these  errors  up  to  ^=4°  is  0.937m2,  while  the  total  rms  including  all  the 
values  in  (4.18)  is  1.869m2. 


Similar  to  the  above,  the  differences  between  the  third  and  first  cases, 
called  the  third  case  errors,  are  obtained  directly  from  (4.9)  as 

(4.19) 


0, 

-9.985, 

-9.291, 

-8.428, 

-7.484, 

-6.507, 

-5.531, 

-4.581, 

-3. 

.675, 

2.824, 

-2.036, 

-1.318, 

-0.674, 

-0.105, 

0.388, 

0.806, 

1.153, 

1. 

.429, 

1.639, 

1.788, 

1.879, 

1.918, 

1.909, 

1.860, 

1.774, 

1.657, 

1. 

.516. 

For  15°  and  30°  these  values  would  be  0.796  and  0.248,  respectively.  The  rms 
of  these  errors  up  to  ^=4 0  is  7.250m2,  while  the  total  rms  including  all  the 
values  in  (4.19)  is  4.240m2.  It  is  thus  apparent  that  the  second  case  is  far 
superior  to  the  third  case  with  regard  to  accuracy,  the  first  case  representing 
the  standard  of  comparison.  In  particular,  the  rms  of  the  second  case  errors 
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attains  only  about  44$  of  the  rms  of  the  third  case  errors.  The  comparison 
of  the  solutions  in  a  specific  exarrple  that  follows  will  confirm  the  superior 
quality  of  the  second  case  (versus  the  third  case)  from  yet  another  angle. 
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4.5  Practical  Example 


In  this  example  only  m*13  geoid  undulations  considered  as  observations 
are  used,  distributed  in  1°  intervals  and  spanning  a  12°-arc.  Upon  taking 
every  other  entry  of  (4.9)  into  account,  the  first  case  is  now  represented 
by 


10.476 

9.291 

7.484  ... 

-1.909 

-1.774 

9.291 

10.476 

9.291  ... 

-1.879 

-1.909 

cl- 

7.484 

9.291 

10.476  ... 

-1.639 

-1.879 

-1.774 

-1.909 

-1.879  ... 

9.291 

10.476 

In  the  second  case,  every  other  entry  of  (4.11)  similarly  entails 


10.476 

8.087 

6.243  ... 

0.608 

0.469 

8.087 

10.476 

8. 087  . . . 

0.787 

0.608 

6.243 

8.087 

10.476  ... 

1.020 

0.787 

0.469 

0.608 

0.787  ... 

8.087 

10.476 
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2  2 

The  new  p  is  now  0.7719  (i.e.,  p  from  4.10b)  and  the  new  p  is  0.5958.  In 


analogy  to  the  procedure  that  yielded  (4.12)  and  (4.13)  one  obtains 


(4.21) 


0.2362 

-0.1823 

0 

0  . . . . 

0 

0 

-0.1823 

0.3769 

-0.1823 

0 

0 

0 

0 

-0.1823 

0.3769 

-0.1823  .... 

0 

0 

0  .  -0.1823  0.2362 


Finally,  the  form  of  the  C  and  P  matrices  in  the  third  case  is 


C3  =  10.476  I  , 


P3  =  0.09546  I  , 


(4.22) 


which  corresponds  to  (4.14)  and  the  equation  that  preceded  it,  except  that 
now  m=13. 

The  13  observations  considered  are  used  to  adjust  one  parameter,  taken 
as  a  point-mass  magnitude.  The  point  mass  itself  is  chosen  to  have  a  central 
location  with  respect  to  the  arc,  and  to  be  located  at  the  depth  d=  1,144km 
underneath  the  earth's  surface  represented  by  a  sphere  of  radius  R=  6,371km. 
Thus  the  central  angle  ^  varies  between  0  and  6°  on  either  sid«=>  ->f  the  point 
mass.  The  mathematical  model  for  geoid  undulations  can  now  be  written  as 


N.  =  (l/G)(kM).M..  ,  i=  1,2,  ...  m  , 


(4.23) 


where  G  is  the  average  value  of  gravity  adopted  as  980  gal;  (kM) ^  =  x  is  the 


! 


unknown  parameter,  i.e.,  the  gravitational  constant  times  the  mass  of  the 
point  mass  "j";  and  «,  is  the  distance  between  the  observational  point  "i" 
and  the  point  mass  "j".  The  latter  is  computed  as 

Aij  =  {[(R-d)sint|>.  ,]2  +  [R-(R-d)cos^Aj  f}1/2  ,  (4.24) 

where  ^i_.  is  the  angle  ip  corresponding  to  the  above  points  i  and  j. 

The  13  observations  of  geoid  undulations  (in  meters)  are  simulated  in  a 
symmetrical  manner,  the  first  and  the  last  being  zero  and  the  middle  one 
directly  above  the  point  mass  being  3.5m,  as  follows: 

0.000,  0.906,  1.750,  2.475,  3.031,  3.381,  3.500, 

3.381,  3.031,  2.475,  1.750,  0.906,  0.000.  (4.25) 

These  values  represent  the  elements  of  the  vector  e  encountered  earlier.  The 
units  are  chosen  such  that  x  is  expressed  in  terms  of  10~6  times  the  earth’s 
kM  (the  latter  is  3.986005  *1014m3/ sec2).  The  matrix  A  has  now  the  dimensions 
(1x13);  the  values  in  this  row,  symmetric  about  its  middle  element,  are 

31.440,  32.542,  33.536,  34.375,  35.015,  35.417,  35.554, 

35.417,  35.015,  34.375,  33.536,  32.542,  31.440. 

(4.26) 

Similar  to  (4.15),  the  least-squares  adjustment  yields 
x  =  (ATPA)_1ATPe  ,  (4.27a) 

o\  =  (AtPA)-1  ,  (4.27b) 

which  entails  three  cases,  P=Plt  P=P2  and  P=P3  ,  in  agreement  with  the  previous 
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classifications.  In  terms  of  the  transformation  x'  =  lOOx  (x1  is  thus  expressed 
in  10-8  of  the  earth's  kM),  the  pertinent  results  can  be  grouped  in  Table  1, 
which  is  self-explanatory.  This  table  shows  the  error  in  the  tri-diagonal 
solution  as  being  only  39%  of  the  error  in  the  diagonal  solution.  The  error 
in  the  sigma  of  the  tri-diagonal  solution  is  even  smaller,  amounting  to  merely 
22%  of  the  error  in  the  sigma  of  the  diagonal  solution.  Although  this  example 
is  quite  specialized,  it  nevertheless  illustrates  that  the  tri-diagonal  set¬ 
up  is  likely  to  match  the  rigorous  setup  much  more  closely  than  would  be 
possible  with  the  diagonal  setup. 


Case 

Name 

A  > 
x' 

Error 

°x* 

Error 

1 

"rigor." 

3.675 

— 

5.491 

— 

2 

tri-diag. 

4.651 

+0.976 

6.121 

+0.630 

3 

diagonal 

6.179 

+2.504 

2.649 

-2.842 

Table  1* 

The  least-squares  solution  and  its  sigma,  together  with  their 
errors,  in  the  three  weighting  schemes  considered;  the  first 
case  is  taken  as  an  errorless  standard. 
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5.  LARGE  AREA  ADJUSTMENT  WITH  POINT-MASS  AND  TIDAL  PARAMETERS 
BASED  ON  THE  BANDED-BORDERED  STRUCTURE  OF  NORMAL  EQUATIONS 

5.1  General  Discussion 

This  chapter  is  a  natural  continuation  of  Chapter  3  in  [Blaha,  1983], 
where  a  large  area  adjustment  algorithm  was  developed  in  terms  of  point-mass 
(P.M. )  magnitudes  as  parameters.  That  algorithm  was  characterized  by  a 
banded  structure  of  normal  equations.  It  was  based  on  a  first-phase  global 
adjustment  of  satellite  altimetry  which,  optionally,  could  have  '.ncluded  also 
certain  tidal  parameters.  In  the  next  section  the  above  algorithm  will  be 
extended,  in  the  sense  that  selected  tidal  parameters  will  optionally  be 
subject  to  (re) adjustment  also  in  the  second-phase  adjustment  of  satellite 
altimetry  based  on  point  masses. 

Since  the  only  modification  of  the  previous  P.M.  algorithm  consists  in 
adding  the  parameters  expressing  certain  temporal  variations  in  the  ocean 
surface,  the  considerations  of  the  underlying  gravity  field  of  the  earth  are 
unchanged  from  the  previous  AFGL  report  ([Blaha,  1983])  to  the  present  one. 

It  may  then  be  useful  to  review  this  approach  as  it  relates  to  the  gravity 
field,  and  compare  it  with  the  P.M.  approach  developed  in  recent  years  and 
described  in  the  paper  [Kautzleben  and  Barthelmes,  1983],  The  former  will  be 
called  "our  approach"  while  the  latter  will  be  abbreviated  as  [KB]  in  the 
forthcoming  comparisons. 

In  our  approach,  the  global  gravity  field  of  the  earth  is  split  into  the 
field  "A",  called  the  global  normal  field,  and  the  field  "B",  called  the 
anomalistic  field.  The  field  A  is  not  an  ellipsoidal  field,  but  a  field 
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represented  by  a  (14,14)  spherical -harmonic  expansion  of  the  earth's  potential; 
all  the  anomalistic  geophysical  quantities,  such  as  geoid  undulations,  gravity 
anomalies,  deflections  of  the  vertical,  etc.,  refer  to  this  field.  The  field 
B  is  described  by  an  equilateral  grid  of  point  masses  at  a  certain  depth 
underneath  the  spherical  surface  approximating  the  earth's  surface  (due  to 
the  relatively  high  resolution  power  of  the  field  A  as  compared  to  the  ellip¬ 
soidal  field,  the  spherical  approximation  is  acceptable). 

In  the  [KB]  approach  the  gravity  field  is  similarly  split  as  can  be 
gathered,  e.g.,  from  the  statement  on  page  308: 

"...  the  anomalies  refer  to  a  global  normal  field  given  by  a  spherical 
harmonic  analysis  up  to  the  fourth  order  inclusively.  Splitting  up  the 
gravity  field  in  a  normal  field  and  the  field  of  anomalies  is  a  mathe¬ 
matically  exact  procedure." 

The  basic  difference  between  the  two  approaches  stems  merely  from  the  fourteenth 
order  spherical -harmonic  expansion  characteristic  of  our  approach  contrasted  to 
the  fourth  order  expansion  of  [KB].  This  difference  is  more  of  a  practical 
than  conceptual  nature. 

Initially,  [KB]  addressed  the  problem  where  not  only  the  P.M.  magnitudes 
were  subject  to  a  least-squares  adjustment,  but  also  their  positions.  Further¬ 
more,  their  number  was  not  fixed.  In  our  approach  the  number  of  point  masses 
as  well  as  their  locations  are  decided  upon  beforehand;  according  to  the 
above  definition  they  are  distributed  in  an  equilateral  grid.  However,  [KB] 
comes  closer  to  our  approach  through  the  text  on  page  310: 

"There  is  a  way  known  [that]  uses  only  points  which  are  distributed 
equidi stantly  on  the  surface  of  a  sphere  in  the  Earth's  interior  and 
where  the  problem  is  reduced  to  determine  the  best  radius  of  that 
sphere  and  [the  magnitudes]  of  the  masses  of  the  points  mentioned." 


-40- 


It  thus  follows  that  in  addition  to  the  solution  of  P.M.  magnitudes  as  in 
our  approach,  [KB]  also  determines  one  additional  parameter,  the  depth  of 
these  masses.  The  common  feature  of  both  approaches  is  the  same  depth  of 
all  the  point  masses  present  in  the  adjustment. 

The  fixed  depth  of  point  masses  in  our  approach  implies  that  the  entire 
adjustment  model  is  linear  in  the  parameters.  This  allows  us  to  adopt  zero 
starting  values  for  the  parameters  --  which  are  unknown  a  priori  —  without 
the  need  for  an  iterated  solution.  On  the  other  hand,  as  it  transpires  from 
[KB]  the  depth  (or  radius)  of  the  point  masses  is  a  non-linear  parameter. 
Therefore,  its  good  starting  value  should  be  known  a  priori  in  order  to  avoid 
iterations  in  the  adjustment  and/or  avoid  the  solution  converging  to  a  false 
minimum.  Clearly,  one  could  perform  preliminary  tests  using  different  values 
of  the  depth  parameter  and  find  the  most  suitable  starting  value  for  the  full 
scale  adjustment.  But  this  can  again  be  related  to  our  approach,  where  dif¬ 
ferent  depths  of  point  masses  can  be  entered  into  adjustment  runs  carried 
out  over  one  and  the  same  limited  area.  The  depth  resulting  in  the  best  fit 
can  then  be  adopted  for  the  desired  full  scale  adjustment.  Accordingly, 
whether  or  not  the  one  additional  parameter,  the  depth  of  the  point  masses, 
should  also  be  subject  to  adjustment  is  essentially  a  practical  matter. 

A  variation  of  the  P.M.  approach  addressed  in  [KB]  is  well  worth  mention¬ 
ing.  It  is  based  on  successive  solutions  of  point  masses,  starting  with  two 
of  them,  as  is  described  on  page  310: 

"After  both  of  these  point  masses  have  been  determined  their  field  is 
subtracted  from  all  the  measured  values.  Based  on  the  new  distribution 
of  the  measured  values  analogously  the  coordinates  and  the  masses  of 
the  next  two  point  masses  are  calculated.  This  algorithm  can  be 
continued  up  to  any  number  of  point  masses." 
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It  is  clear  that  whether  two  or  more  than  two  point  masses  are  treated  at  a 
time  in  this  manner  is  an  arbitrary  decision.  One  could  as  well  consider  all 
of  the  point  masses  and  subtract  their  field  from  the  measured  values  in  one 
step.  This  is,  in  fact,  the  principle  used  in  our  approach,  i.e.,  a 
simultaneous  least-squares  adjustment  of  all  the  P.M.  parameters. 

Our  simultaneous  approach  could  have,  indeed,  some  practical  drawbacks 
from  the  computer-core  and  run-time  standpoints.  However,  the  development  of 
a  banded  (or  banded-bordered)  structure  of  normal  equations  can  bridge  these 
difficulties.  Such  a  structure  results  in  a  ten-fold  increase  in  efficiency, 
and  has  allowed  the  increase  in  the  number  of  P.M.  parameters  adjusted 
simultaneously  from  about  200  to  2,000. 

One  of  the  three  basic  steps  leading  to  the  banded  structure  of  normal 
equations  described  in  Chapter  3  of  [Blaha,  1983]  resides  in  neglecting  a 
point  mass  when  forming  the  observation  equation  for  a  given  observation  if 
this  mass  is  farther  from  the  observation  point  than  a  predetermined  limit. 
But  this  is  compatible  with  the  procedure  used  in  [KB]  on  page  310: 

"To  determine  the  improved  values  for  the  coordinates  and  for  the 
masses  of  both  these  point  masses  by  the  procedure  described  above 
not  all  the  measured  values  at  the  whole  surface  of  the  Earth  but  only 
values  in  the  neighbourhood  of  these  point  masses  are  used." 

The  approximation  of  neglecting  an  observation  when  its  influence  on  a  point 
mass  is  only  a  fraction  of  the  maximum  influence  (for  an  observation  directly 
above  the  point  mass)  is  well  worth  making,  both  in  [KB]  and  in  our  approach. 
One  can  thus  conclude  that  although  some  of  their  individual  features  may 
differ,  the  two  approaches  on  the  whole  are  conceptually  quite  similar. 


-42- 


5.2  Banded-Bordered  Algorithm 


The  task  addressed  in  this  section  is  the  extension  of  the  concept 
and  the  algorithm  for  large-scale  adjustments  of  satellite  altimetry  using 
point  masses.  Whereas  in  Chapter  3  of  [Blaha,  1983]  special  attention  was 
paid  to  the  P.M.  algorithm  applicable  to  a  banded  system  of  normal  equations 
eventually  resolved  by  the  "modified  Choleski  algorithm",  the  present 
extension  concerns  a  banded-bordered  system  of  normal  equations.  It  will 
become  apparent  that  the  former  system  is  a  special  case  of  the  latter. 

First,  the  highlights  of  the  resolution  of  the  banded  system  are  re¬ 
capitulated.  The  observation  equations  in  usual  notations  read 

V  -  AX  -  (Lb-L°)  ,  (5.1) 

where 

A  =  (design)  matrix  of  observation  equations, 

X  =  vector  of  parameters  (here  the  P.M.  magnitudes), 

_(Lb-L°)  =  vector  of  constant  terms  of  observation  equations, 

and  where 

Lb  =  observed  values  of  the  observables, 

L°  =  initial  values  of  the  observables. 


The  normal  equations  for  the  above  system  are 

NX  =  U  ,  (5.2) 

where 

N  =  AtA  ,  U  =  AT(Lb-L°)  .  (5.3a,b) 

The  weight  matrix  of  the  observations,,  P,  is  not  considered  here  because  the 
matrix  A  is  scaled;  at  this  point  one  would  simply  have  P=I  (a  unit  matrix). 

In  the  case  of  independent  observations  the  scaling  is  accomplished  by 
dividing  each  observation  equation  by  the  appropriate  sigma,  which  was  used 
in  [Blaha,  1983  ]  and  mentioned  on  page  33  therein.  In  the  situation  with 
weighted  parameters  the  appropriate  weight  matrix  Px  would  be  added  to  the 
matrix  of  normal  equations,  N.  However,  the  initial  values  of  the  P.M.  para¬ 
meters  (taken  as  zeros)  are  not  attributed  any  weight. 

The  solution  of  normal  equations  and  the  variance-covariance  matrix  of 
the  parameters  are  expressed  by 

X  =  N-1U  ,  E  =  N-1  .  (5.4a,b) 

A 

With  regard  to  a  w-vector  of  linear  combinations  of  parameters,  we  have 

F  =  U^X  ,  ,  (5.5a,b) 

where  U  is  a  matrix  of  dimensions  (n*w),  n  being  the  number  of  parameters. 

The  variance-covariance  matrix  of  the  parameters  is  obtained  by  substituting 
I  for  U  in  (5.5b),  which  is  the  procedure  adopted  in  the  modified  Choleski 
algorithm. 
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As  is  explained  on  page  37  of  [Blaha,  1983],  the  triangular  decomposition 
of  the  positive-definite  matrix  N  amounts  to  writing 


TTT  =  N  , 

N'1  =  TV)-1  ; 

(5.6a,b) 

ttr  =  U  , 

73 

II 

1 

<= 

(5.7a,b) 

T1^  =  U  , 

R  =  (tt)"1u  . 

(5.8a,b) 

The  matrix  T  of  dimensions  (n*n)  is  computed  according  to  the  explanation 
given  on  pages  38  and  39  of  the  above  reference,  listing  also  the  remaining 
part  of  the  algorithm  ("regular"  Choleski  algorithm  at  this  point)  for  the 
elements  of  R,  R  and  X.  The  latter  is  computed  via  the  "back  substitution", 
corresponding  to 

TX  =  R  ,  X  =  T-1R  ,  (5.9a,b) 

where  (5.9a)  is  the  consequence  of  (5.2)  being  pre-multiplied  by  (Tt)-1  , 
and  (5.9b)  follows  either  from  (5.9a)  or  from  (5.4a)  with  (5.6b),  (5.7a). 
Finally,  (5.5b)  can  be  expressed  as 

(5.10) 

which  leads  to  £x  upon  substituting  I  for  U. 

In  the  present  task,  the  vector  A'X'  is  added  to  the  right-hand  side  of 
(5.1),  subject  to  the  same  least-squares  process.  Here  X'  represents  the 
"overcorrections"  to  the  tidal  parameters  which  were  involved  in  the  first- 
phase  global  altimeter  adjustment  together  with  the  spherical -harmonic  poten¬ 
tial  coefficients  and  the  state  vector  parameters;  the  matrix  A'  of  the 
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partial  derivatives  with  respect  to  the  tidal  parameters  is  adopted  from 
that  adjustment  (its  elements  are  stored  on  a  tape  for  each  altimeter 
observation).  The  initial  values  for  X'  are  set  to  zero  and  the  weight 
matrix  for  X',  denoted  P  ,  reflects  the  reliability  of  the  previous 
solution  for  the  tidal  parameters,  i.e.,  the  strength  of  the  zero  initial 
values.  The  extended  observation  equations  read 

V  =  [A  A'] fin  -  (Lb-L°)  .  (5.11) 

X* 

—  — 


The  extended  normal  equations  then  are 


with  the  new  notations. 


U  =  AtA'  ,  N'  «  A,tA'  +  P  ,  U'  =  A,T(Lb-L°)  . 

A 


(5.12) 


(5. 13a, b,c) 


The  weight  matrix  P  is  incorporated  in  A'  through  the  same  scaling  as  before. 

The  solution  of  the  above  normal  equations  and  the  variance-covariance 
matrix  of  the  parameters  are  given,  respectively,  as 

N  t? 

IF  N' 


Upon  using  equation  (11a)  from  [Blaha,  1976],  one  has 


— 

__  I 

-1 

-  _  _  -V 

N 

N-1  +  n_1uhFn_1  -n_1uh 

F 

a* 

N' 

mg 

-hFn-1  h 

where 

H  =  (N'  -  FtT1!)-1  .  (5.16) 

With  (5.15)  and  (5.16),  equations  (5.14a,b)  lead  to 

X'  =  (N‘  -  UtN_1U)‘1(U'  -FlCHj)  ,  (5.17a) 

X  «  N-1(U  -  lx')  ;  (5.17b) 

Ex,  «  (N‘  -  Fn-1^)"1  ,  (5.18a) 

Ev  =  N-1  +  N-1U(N'  -  UTN“1TJ)-1UTN“1  .  (5.18b) 

A 


In  analogy  to  (5.5a,b),  we  write  for  linear  functions  of  X'  and  X, 
respectively: 

F'  =Tj,TX'  ,  IF.  =FTZx,Ft  (5. 19a, b) 

F  =  Fx  ,  =  tFyJ  ,  (5.20a,b) 

where  the  vector  F'  has  w'  elements  and  TJ'  is  a  matrix  of  dimensions  (n'xw'), 
n'  being  the  number  of  tidal  parameters  (i.e.,  the  number  of  elements  in  X'). 


F 


The  variance-covariance  matrices  forX'  and  X  are  obtained  upon  setting 
U 1 = I  in  (5.19b)  and  U=I  in  (5.20b),  respectively. 

As  an  extension  of  (5.7a,b)  and  (5.8a,b),  new  notations  are  introduced, 
namely 

T1!  =  IT  ,  =  (TT)_1f  .  (5.21a,b) 

Thus  the  following  equivalent  expressions  are  obtained: 

OV1!!  =  Ff  ,  ^>rN_1U  =  Fr  ,  N_1f  =  T"1!  . 


Accordingly,  (5. 17a, b) , (5. 19b)  and  (5.20b)  are  rewritten  as 
X'  =  (N1  -  Fr)"1^'  -  Fr)  , 

X  =  T—1 ( R  -  ffX')  ; 

Ef,=  IPT(N'  -  Fr")'1!?  , 

E  =  Fr  +  F^N1  -  Ff)-1^  , 

F 


(5.22a) 

(5.22b) 

(5.23a) 

(5.23b) 


where,  in  addition  to  (5.21b),  also  the  previous  notations  (5.6b)  and  (5.8b) 
have  been  used. 


Since  H  in  (5.15)  and  thus  also  H"1=N’-FR:  are  positive-definite  matrices, 
one  can  express  X'  in  (5.22a)  and  E_,  in  (5.23a)  in  a  complete  analogy  to  the 
process  represented  by  the  steps  (5.6a)-(5. 10).  In  particular,  we  write 


T,tT' 

=  N’ 

-  pr , 

(N’  -  Ft)'1  =  T'"1(T't)"1  ; 

(5.24a,b) 

T,TR' 

=  U' 

-  Pr  , 

R*  =  (rW  -  Fr)  ; 

(5.25a,b) 

T'1!’ 

=  U’ 

» 

I1  «  (t,t)"i:u'  . 

(5.26a,b) 

i 


t 
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This  then  leads  to 


T*X*  =  R'  ,  X1  =  T,_1R'  ,  (5.27a, b) 

IF,  -  f'T^'  .  (5.28) 

Upon  adding  one  further  notation,  R',  and  using  it  in  the  relation 

— rp  _  m  ^  1 — /p 

T'  R'  =  R  R  ,  R'  =  (T‘  )  R  R  ,  (5.29a,b) 

it  follows  from  (5.23b)  that 

Ef  =  RTR  +  R,TR'  .  (5.30) 

The  whole  process  can  now  be  recapitulated  in  six  steps  as  follows. 

1)  Equations  (5.6a),  (5.7a),  (5.8a)  and  (5.21a)  describe  the  part  of 
the  algorithm  that  can  be  schematically  represented  as 

TT[T|RjR;I]=[N;UjUjU], 
nxn  i  nxl  i  nxw  •  nxiV  nxn  •  nxi !  n*w  I  nxn^ 

where  the  dimensions  are  also  indicated  and  where  the  matrices  under¬ 
lined  by  a  dashed  line  pertain  to  the  new,  "border"  part  of  the 
algorithm.  From  here  T,  R,  R  and  R  are  all  computed  as  in  the  process 
mentioned  following  (5.8a,b). 

2)  Equations  (5.24a),  (5.25a),  (5.29a)  and  (5.26a),  respectively, 
describe  the  part 

T,t[  T'  1  R'  j  R'  {  R'  ]  -  C(N'-Fr)  !  (u'-Fr)  [Fr  j  U'  ]  , 

n'xn"  n'xi  'n’xw  'n'xw'  n'xn'  1  n'xi  'n'xw'n'xw' 
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which  pertains  to  the  "border"  part  of  the  algorithm  in  its 
entirety.  From  here  T’,  R1,  ft'  and  R-'  are  computed  in  the  same 
way  as  their  unprimed  counterparts  above. 

3)  Equation  (5.27a)  is  repeated  as 

T'  X’  =  R’  , 

(VxjV  __n]_xl _ n^i 

which  again  pertains  to  the  "border"  part  of  the  algorithm.  It 
yields  X'  upon  the  "back  substitution". 

4)  Equation  (5.22b)  is  rewritten  as 

T  X  »  (  R  -  R  X')  , 

n*n  nxl  nxl  nxl 


from  which  X  is  obtained  by  the  "back  substitution". 

5)  Equation  (5.30)  is  rewritten  as 
Ef  =  +  tf,TE'  . 

WXW  WXW  WXW 


6)  Finally,  equation  (5.28)  is  rewritten  as 

zpI  =  . 

W^*W'  w'xw' 

All  of  the  steps  above  must  be  carried  out  in  the  order  indicated,  except 
perhaps  for  the  steps  5)  and  6).  The  banded  system  is  seen  as  a  special  case 
of  the  banded -bordered  system  above,  the  additions  due  to  the  "border"  being 
identified  by  the  dashed  lines.  Thus  the  banded  system  as  described  in  [Blaha, 
1983]  consists  only  of  the  steps  1),  4)  and  5)  with  the  self-evident  sinpllfica- 
tions.  Advantage  of  the  banded  structure  of  normal  equations  —  in  addition  to 


the  advantage  of  N  being  positive-definite  —  can  be  taken  in  the  step  1)  as 
described  in  the  above  reference,  where  N  and  T  are  compacted  to  form  (bxn) 
matrices,  b  being  the  bandwidth;  this  reference  uses  the  symbol  T  for  the 
modified  matrix  N  while  the  symbol  T  is  left  unchanged. 

The  steps  1 ) -4 )  above  are  now  depicted  symbolically  with  regard  to  the 
replacements  in  the  computer  core  as  they  are  progressing  in  the  course  of  the 
present  "extended  modified  Choleski  algorithm": 

CT]CU!U!U]-CT][R!RJR], 

bxn  nxi  ■  nxw  1  nxn1  bxn  nxi  •  nxw  1  nxn' 


[N’-R‘rR!U'-RTR,!RTR  J  U' 
n'xn'  n'xl  n'xw  n'xw' 


]  +  [T'  •  R'  i  R'  j  R’  ] 
n'xn'  n'xl  n'xw  n'xw' 


r  1  I  __  I  “Z 

-*■  [  T'  !  X'  1  I'  J  If'  ], 
n'xn'  n'xl  n'xw  n'xw* 


[  T  ][  R  -  RX'  R  5  R  ]  -*■  [  T  ][  X  i  R  'J  R  ]  . 

bxn  nxi  nxi  nxw  nxn^  bxn  nxi  nxw  nxn' 

The  steps  5)  and  6)  need  not  be  listed  again.  This  scheme  provides  for  the 
resolution  of  the  banded-bordered  systems  of  normal  equations.  Here  again, 
if  the  underlined  parts  are  left  out,  one  is  in  the  presence  of  the  algorithm 
for  the  banded  systems  as  described  in  Chapter  3  of  [Blaha,  1983]. 
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6.  LEAST-SQUARES  COLLOCATION  WITH  NOISE  AS  A  SECOND-PHASE 

ALTIMETRIC  APPROACH 

In  Chapter  5  as  well  as  in  [Blaha,  1983]  a  second-phase  geoidal 
determination  has  been  addressed  via  a  "direct"  approach,  in  the  sense  that 
the  observations  subjected  to  the  least-squares  (L.S.)  adjustment  have 
been  treated  at  their  original  locations.  The  quotes  are  used  to  indicate 
a  special  point  of  view,  since  in  another  aspect  such  an  adjustment  could 
be  considered  as  indirect.  In  particular,  it  proceeds  to  an  evaluation  of 
the  desired  quantities  describing  the  earth's  gravity  field  (geoid  undulations, 
gravity  anomalies,  etc.)  through  intermediate  quantities,  namely  the  point 
masses.  Due  to  the  simultaneous  character  of  this  adjustment  with  regard 
to  both  the  observations  (minus  the  geoidal  residuals  at  the  original 
locations)  and  the  parameters  (point-mass  magnitudes),  it  is  computationally 
demanding.  This  is  true  especially  for  a  large-scale  adjustment,  in  which 
case  the  area  is  subdivided  into  overlapping  strips  treated  separately. 

In  attempts  for  an  increasingly  detailed  representation  of  the  earth's 
gravity  field  on  a  global  scale  it  thus  becomes  expedient  to  search  for  ways 
which  would  complement  the  "direct"  approach  by  using  other  methods  than  the 
least-squares.  For  example,  one  could  take  advantage  of  integral  formulas 
and  compute  spherical -harmonic  (S.H.)  potential  cofficients  consistent  with 
the  desired  smoothing  effect  without  proceeding  to  a  simultaneous  L.S. 
adjustment  of  observations.  The  final  result  would  be  a  uniform  representa¬ 
tion  of  the  geoid  and  the  earth's  gravity  field  obtained  for  the  whole  globe 
in  a  continuous,  one-piece  fashion.  The  trade-off  would  be  the  necessity 
of  using  some  kind  of  a  "translocation"  approach. 
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Various  interpolation  and  approximation  techniques  exist  where,  unlike 
in  the  point-mass  approach,  the  original  observations  are  translocated  to 
some  strategically  advantageous  positions.  When  used  in  the  role  of  inter¬ 
polating  the  (residual)  geoid  undulations  of  satellite  altimetry,  the  L.S. 
collocation  with  noise  pioneered  by  Moritz  [1980]  could  be  considered  as 
belonging  to  a  "translocation"  category.  A  global  approach  based  on  Moritz's 
collocation  will  -onstitute  the  main  topic  of  this  chapter.  Besides  serving 
in  its  own  right,  the  collocation  approach  could  be  further  exploited  and 
result  in  a  set  of  S.H.  potential  coefficients  as  suggested  in  the  preceding 
paragraph. 

The  basic  observed  quantities  to  be  used  in  the  collocation  approach  are 
described  in  [Blaha,  1983].  They  consist  of  the  "observed"  geoid  undulations 
or,  equivalently,  of  the  corresponding  high-resolution  residuals  which  en¬ 
compass  the  improvements  due  to  the  orbital  and  tidal  adjustment  of  satellite 
altimetry.  Just  as  the  "direct"  approach,  the  collocation  approach  and, 
eventually,  the  S.H.  approach  can  also  serve  in  predicting  other  geophysical 
quantities  in  addition  to  geoid  undulations,  such  as  gravity  anomalies, 
deflections  of  the  vertical,  etc.  Clearly,  the  S.H.  approach  would  have  to 
be  carried  out  on  a  global  scale.  In  all  cases  the  final  product  is  repre¬ 
sented  by  geoidal  contour  maps  and,  optionally,  by  contour  maps  of  other 
geophysical  quantities,  especially  the  gravity  anomalies.  And  in  all  cases, 
the  contour  maps  are  desired  on  a  global  oceanic  scale. 
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6. 1  Simple  Translocation  Approach 

This  section  is  conceived  as  an  illustration  of  a  translocation  approach 
which  is  the  simplest  but  also  the  least  accurate.  In  this  approach,  the 
location  of  an  observation  point  is  changed  but  the  value  of  the  "observed" 
geoid  undulation  is  not  modified  in  any  way.  A  grid  is  generated  by  filling 
it  with  undulation  values,  each  adopted  simply  as  the  value  at  the  nearest 
observation  point.  The  formation  of  a  regional  or  global  grid  of  geoid 
undulations  is  thus  accomplished  very  fast,  from  which  a  contour  map  can  be 
constructed  by  simple  means  such  as  polynomial  interpolation  (linear  or 
higher-degree).  In  considering  the  distribution  of  SEASAT  passes,  one  could 
generate  an  equilateral  grid  of  dimensions  l°xl°  or  larger.  In  order  to 
facilitate  the  interpolation  procedure  in  some  practical  cases,  a  geo¬ 
graphical  grid  may  be  generated  instead. 

The  disadvantages  of  this  approach  are  listed  as 

a)  The  values  associated  with  the  "observed"  geoid  could  be  grossly 
in  error  at  the  grid  points. 

b)  Such  errors  are  unpredictable,  neither  minimized  nor  estimable, 
due  to  the  grid  values  obtained  without  any  regard  to  statistical 
methodology. 

c)  The  contours  may  contain  additional  errors  due  to  the  inter¬ 
polation  procedure  (a  polynomial  interpolation  makes  no  provision 
for  the  physical  characteristics  of  the  quantity  being  inter¬ 
polated);  in  particular,  the  grid  cannot  be  densified  through  a 
physically  meaningful  model  before  the  application  of  the 
polynomial  interpolation. 
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d)  Contour  maps  of  other  geophysically  meaningful  quantities 
(e.g.  gravity  anomalies)  cannot  be  constructed. 

e)  The  geographical  grid,  if  used,  leads  to  an  apparent  increase 
in  resolution  at  higher  latitudes  and  in  the  E-W  direction, 
which  is  completely  artificial  and  in  contradiction  with  the 
physical  reality  (the  geoidal  detail  is  independent  of  latitude 
or  azimuth). 

f)  A  great  majority  of  observations  are  ignored. 


6.2  Review  of  the  Best  Linear  Prediction  Algorithm 


The  approach  just  discussed  can  serve  only  for  comparison  purposes  and 
does  not  represent  the  final  outcome  of  satellite  altimetry  in  any  way.  A 
meaningful  outcome  based  on  altimeter  measurements  can  be  obtained  with  the 
aid  of  the  L.S.  collocation  with  noise,  to  be  described  in  the  next  sections. 

It  will  then  become  apparent  that  the  shortcomings  listed  above  as  a)  through 
f)  disappear.  However,  before  presenting  this  method  we  review  its  special 
case,  called  here  the  best  linear  prediction  method,  whose  detailed  derivation 
can  be  found  in  Appendix  2.  In  [Moritz,  1980]  this  case  is  called  the  error¬ 
less  collocation. 

In  the  context  of  satellite  altimetry  only  one  prediction  may  be  performed 
at  a  time.  Although  the  next  section  will  present  the  algorithm  for  the  L.S. 
collocation  with  noise  in  the  matrix  form,  a  "transliteration"  similar  to  the 
one  below  would  be  exceedingly  simple  to  carry  out.  In  denoting  the  prediction 
point  by  J,  and  the  observation  points  (where  the  geoid  undulations  are  available) 
by  1,  2,  ....  n,  the  following  formulas  in  the  best  linear  prediction  approach 
can  be  transcribed  from  Appendix  2: 

Nj  =  MtW  ,  (6.1) 

where 

M  =  CDji  °j2  “•  °jn3  ’  (6-2) 

w  =  H_1F  ,  (6.3) 
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with 


F  =0i  N2  ...  Nj  .  (6.5) 


In  these  formulas  DQ  =  D(0),  Di^  =  D(4|i ^ )  are  the  values  of  the  covariance 
function  for  geoid  undulations,  N1,  N2,  ....  Nn  are  the  geoid  undulations  at 
the  observation  points,  and  Nj  is  the  desired  prediction  of  geoid  undulation. 
The  dispersion  measure  (vN)  due  exclusively  to  the  error  of  prediction  is  com¬ 
puted  in  this  case  as 


and  it  can  be  viewed  in  a  loose  analogy  to  the  variance  in  stochastic  processes; 
in  Appendix  2  it  is  a  diagonal  element  of  the  dispersion  matrix  M{eeT}  . 

If  gravity  anomalies  are  the  predicted  quantities  (based  on  geoid  undula¬ 
tions  at  the  observation  points),  one  can  write  formulas  similar  to  those  above: 

Aga  =  MTW  ,  (6.6) 

where 

•••  BJ„3-  (6.7) 
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and  where  W  is  the  same  vector  as  before,  while  H_  h  H^..)  is  the  cross¬ 
covariance  function  as  presented  in  Appendix  2.  The  dispersion  measure  is 
now 

v?  =  C  -  , 

Ag  o 

which  corresponds  to  a  diagonal  element  of  M{eeT} ;  Cq  h  C(0)  in  this  expression 
is  obtained  from  the  covariance  function  for  gravity  anomalies. 


6.3  Collocation  Algorithm 


The  formulas  for  the  L.S.  collocation  with  noise  can  be  developed  in  a 
close  analogy  to  the  best  linear  prediction  method  of  Appendix  2.  The  basic 
difference  between  the  two  approaches  resides  in  the  statistical  nature 
attributed  to  the  observations.  In  particular,  in  the  present  approach  the 
observations  are  no  longer  considered  as  errorless  quantities.  Thus,  in 
addition  to  the  symbol  M  (averaging  operator  over  the  unit  sphere)  introduced 
in  the  beginning  paragraph  of  Appendix  2,  also  the  symbol  E  (mathematical 
expectation  operator)  would  have  to  be  considered  in  the  derivations.  In 
adapting  the  development  of  Appendix  2  to  the  present  context,  one  can  imagine 
the  operator  M  replaced  by  the  operator  M.  The  latter  represents  both  of  the 
operators  E  and  M  taken  together  in  either  order  (due  to  their  independence), 
symbolized  by 

M  =  EM  =  ME  . 

With  the  above  stipulation  the  formulas  of  Appendix  2  can  be  easily 
modified.  In  particular,  the  errorless  vector  F  in  (A2.7-9)  is  replaced  by 
Fb,  where 

Fb  =  F  +  e  , 

and  where  e  is  the  vector  of  observational  noise  whose  variance-covariance  matrix 
is 

E(eeT}  s  e  . 
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From  the  nature  of  the  operators  E  and  M,  described  in  [Moritz,  1980],  it  follows 
that 


E (P }  *  P  ,  E{F}  *  F  ,  E{e)  =  0  , 

M{P>  =  0  ,  M{F}  =  0  ,  M{e}  =  e  . 

As  further  consequences,  we  have 

M{PPT)  =  M{PPt}  , 

M(FbPT}  =  M{FPt)  +  M{ePT) 

5  M{FPt}  +  E{e}M{PT}  =  M{FPt}  , 

M{fbFbT}  =  M{FFt}  +  M(FeT)  +  M{eFT}  +  M{eeT) 

=  M{FFt}  +  M{eeT}  =  M{FFt}  +  E(eeT}  .  (6.8) 

With  the  above  changes,  (AZ.9)  in  the  present  context  would  read  as  it 
stands  except  that  M  on  the  left-hand  side  would  be  replaced  by  M,  and  M{FFT} 
on  the  right-hand  side  would  have  E{eeT)  added  to  it.  This  last  change  would 
take  place  also  in  (A2. 10-12);  in  addition,  the  vector  F  in  (A2.ll)  would  be 
replaced  by  Fb  and  the  operator  M  on  the  left-hand  side  of  (A2.12)  would  be 
replaced  by  M.  Upon  denoting  the  "error  measure"  as  described  through  M  just 
mentioned  by  the  symbol  C,  the  final  results  corresponding  to  (A2.ll1),  (A2.12') 
and  the  remainder  of  that  paragraph  in  Appendix  2  would  now  read 


P  =  MT(H+  E)_1Fb  =  MTW  ,  (6.9) 

C.  =  Hp  -  MT(H+E)_1M  ,  (6.10) 
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where 


W  =  (H+  E)_1Fb  ; 

H  =  M{FFt}  ,  E  e  E(eeT}  ; 

Mt=  M{PFt}  ,  HpE  M{PPt}  . 

Equations  (6.9)  and  (6.10)  are  essentially  the  formulas  (14-27)  and  (14-42)  in 
[Moritz,  1980],  characterizing  the  L.S.  collocation  with  noise. 

With  regard  to  the  predictions  P  symbolizing  functions  of  different  kind 
from  F,  the  results  corresponding  to  (A2.11")  and  (A2.12")  in  Appendix  2  would 
similarly  become 

A 

P  =  MT(H  +  E  r1^  =  MTW  ,  (6.11) 

C~  =  Hg  -  M^H  +  E )-1M  ,  (6.12) 

with  the  additional  notations 

MT  e  M{PFt}  ,  =  MCPPT>  . 

The  discussion  in  the  paragraphs  following  (A2.15)  in  Appendix  2  would 

now  be  somewhat  changed.  If  all  the  prediction  points  coincided  with  the  (n) 

observation  points,  it  would  similarly  hold  true  that 

if  =  H  ,  Hp  =  H  ; 

but  from  (6.9)  and  (6.10),  upon  changing  "P"  to  it  would  now  follow: 
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F  =  H(H  +  l)_1Fb  , 


(6.13) 


C^=  H  -  H(H+  E)-1H  .  (6.14) 

u 

These  equations  no  longer  yield  F  and  0*  respectively,  as  in  Appendix  2,  but 
the  results  will  approach  such  values  as  z  approaches  0. 

In  the  case  where  one  prediction  point  (J)  coincides  with  the  one  and 
only  observation  point  (when  n=l),  equations  (6.13)  and  (6.14)  become 

F,  ■  tDo/(Do  +  02)lF^  , 

4  '  D<,-Do/(Doto2)  > 

«J 

where  the  variance-covariance  matrix  E  reduces  to  the  variance  a2.  If  o2«D  , 

O 

the  above  relations  reveal  that 


J 


Thus  F  is  not  reproduced  exactly,  but  is  reduced  by  the  (small)  factor  a  / D  . 

J  O 

Here  again,  as  a 2  goes  to  zero  the  previous  results  are  recovered. 

Instead  of  being  based  on  the  S.H.  potential  coefficients,  the  formulas 
for  cn  ,  dn  and  hn  in  Appendix  2  can  be  developed  with  the  aid  of  the  closed- 
form  expressions  from  [Tscherning  and  Rapp,  1974],  adapted  on  pages  76  and  77  of 
[Blaha,  1982],  For  the  k-th  degree  variances  we  then  have 

ck  =  sk+2A(k-l)/[(k-2)(k+B)]  , 

dk  =  ckR2/[G2(k-l)2]  , 

hk  =  ckR/[G(k-l)]  =  dkG(k-l)/R  , 
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s  =  0.999617  , 


A  =  425.28  mgal 2  , 

B  =  24  , 

and  with  R  =  6.371xl06m  representing  the  earth's  mean  radius  and  G  =  979.8  gal 
representing  the  average  value  of  gravity  on  the  earth's  surface. 

In  the  final  formulas  below  the  degree  n  represents  the  smallest  features 
of  the  earth's  gravity  field  taken  into  account  by  the  known  S.H.  expansion  and 
C,  D,  H  represent  the  (cross-)  covariance  functions  due  to  the  neglected  degrees: 

ck  =  0.999617k+2x(k-l)x425.28  mgal2/ [(k-2)(k+24)]  ,  (6.15) 

dk  =  0.999617k+2xl7,981  m2/ [( k- 1 ) ( k-2 ) ( k+24 ) ]  ,  (6.16) 

hk  =  0. 999617k+2x2 ,765.3  mgal xm/[(k-2) (k+24) ]  ;  (6.17) 


C(*)  =  l 

k=n+l 

ckPk(cos^)  , 

C(0)  « 

00 

l  ck  , 

k=n+l  K 

(6.18) 

D<*)  *  I 

k=n+l 

dkPk(cos^)  , 

D(0)  = 

00 

l  \  ‘ 

k=n+l 

(6.19) 

H(*)  -  l 

k=n+l 

hkpk(c°s^)  , 

H(  0)  = 

00 

I  * 

k=n+l  k 

(6.20) 

We  note  that  (6.18-20)  could  be  computed  using,  for  example,  the  closed 
forms  appearing  on  pages  44,  45  of  [Tscherning  and  Rapp,  1974].  But  if  the 
practical  computation  proceeds  as  indicated  in  (6.18-20),  the  symbol  "<*>"  can  be 
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replaced  by  a  sufficiently  large  number  such  as  500  or  1,000,  etc.  With  the 
above  values  of  the  functions  C,D  and  H,  the  matrices  H,  MT  and  Hp  needed  in 
(6.9)  and  (6.10),  as  well  as  the  matrices  MT  and  H~  needed  in  (6.11)  and  (6.12), 
can  be  formed  as  in  Appendix  2,  equations  (A2. 13-15)  including  the  remainder  of 
that  paragraph. 

The  theory  just  described  is  first  considered  in  view  of  the  predictions 
of  geoid  undulations,  where  the  observed  quantities,  Fb,  are  the  geoid  undulations 
referring  to  the  S.H.  model  truncated  at  (n,n)  degree  and  order.  These  observa¬ 
tions  are  typically  the  (minus)  geoidal  residuals  from  a  S.H.  adjustment  of 
satellite  altimetry.  The  formulas  (6.9)  and  (6.10)  are  now  rewritten  as 

P  =  MTH-1Fb  , 

Cg  =  Hp  -  , 

where 

Fb  =  F+e  ,  (6.23) 

H  =  H+  E  .  (6.24) 

Here  F  is  the  vector  of  "true"  values  of  geoid  undulations  referring  to  the 
(n,n)  field,  i.e.,  of  values  owing  their  existence  to  the  truncated  degrees  from 
n+1  to  infinity,  while  e  is  the  observational  noise  (independent  of  F).  The 
matrix  H  corresponds  to  the  vector  F,  in  the  sense  that  "n"  in  the  covariance 
function  D(^)  used  in  forming  H  according  to  the  locations  of  the  observation 
points  is  the  same  as  the  "n"  of  the  (n,n)  field  associated  with  the  values  In 
F.  In  the  same  vein  the  matrix  MT  corresponds  to  F  as  well,  while  the  matrix 


(6.21) 

(6.22) 
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I  corresponds  to  the  vector  e.  The  matrix  Hp  corresponds  to  P  and  thus  to 
the  locations  of  the  prediction  points,  but  D(^)  and  "n"  used  in  its  con¬ 
struction  are  the  same  as  above. 

Next  consider  a  practical  case  in  which  the  predictions  are  required 
to  describe  a  smooth  geoid.  Such  a  geoid  can  be  represented  by  a  S.H.  expansion 
(n',n'),  as  opposed  to  the  theoretical  expansion  («,«>).  For  example,  a  1°- 
geoid  would  correspond  roughly  to  a  (180,180)  S.H.  expansion.  In  this  task, 
the  "true"  geoid  is  assumed  to  be  described  perfectly  by  a  (n '  ,n ' )  S.H. 
expansion,  all  the  other  S.H.  potential  coefficients  being  considered  as  zeros. 
With  regard  to  such  a  geoid,  the  above  formulas  would  be  unchanged  except  that 
the  symbol  «>  in  the  expression  for  D(i/>)  would  be  replaced  by  n'.  However,  an 
inconsistency  would  manifest  itself  in  a  real-world  situation  where  F  refers  to 
the  actual  (unsmooth)  geoid  as  sensed  through  the  observations  Fb.  Therefore, 
the  adaptation  of  D(^)  to  the  "true"  geoid  under  consideration  is  equivalent  to 
modeling  only  the  corresponding  F'-part  of  F  as  errorless  values  and  pushing 
the  remaining  part,  F",  into  the  realm  of  "errors",  namely 

F  =  F’  +  F"  .  (6.25) 

The  observational  vector  Is  thus  split  as  follows: 

Fb  =  F'  +  (F"  +  e)  ,  (6.26) 

or 

Fb  =  F*  +  e'  ,  (6.27) 

e'  =  F" +  e  .  (6.28) 
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The  problem  of  the  L.S.  collocation  with  noise  is  now  reformulated  in  a 
slightly  modified  fashion.  One  is  to  find  the  best  linear  estimate  P'  , 

P'  =  A'TFb  ,  (6.29) 

where  the  matrix  A'  is  unknown.  The  primes  indicate  that  the  "true"  geoid  cor¬ 
responds  to  the  highest  degree  n'  and  not  to  °°  which  would  characterize  the 
unprimed  quantities.  In  analogy  to  Appendix  2,  the  problem  of  finding  A,T  is 

_  A 

addressed  through  the  operator  M.  With  e'  = P'-P‘,  one  forms  the  "dispersion 
matrix" , 

e'e,T  =  P'P,T  -  P * FbTA 1  -  A,TFbP,T  +  A’W^A'  , 
and  its  "average"  as 

M{e'e,T}  =  M{P,P't}  -  MtP’F^M*  -  A'TM{FbP,T}  +  A,TM{FbFbT}A'  .  (6.30) 

Upon  considering  the  nature  of  the  operators  M  and  E  as  in  [Moritz,  1980], 
it  follows  for  the  first  term  on  the  right-hand  side  of  (6.30): 

M{P'P,t}  =  M{P'P,t}  .  (6.31) 

The  second  term  (and  thus  also  the  third  term)  is  first  developed  according  to 
(6.27)  as 

M{P'FbT)  =  M{P'F,t}  +  M{P'e'T}  .  (6.32) 

Similar  to  (6.31),  we  have 

M{P'F,t}  =  M{P'F,t}  .  (6.33a) 
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On  the  other  hand,  in  considering  (6.28)  it  follows  that 


M(P'e'T}  =  M{P'F"t}  +  M{P‘eT}  =  M{P'F',T}  +  M(P'}E{eT}  h  M{P'F"T}  . 

But  this  last  matrix  is  zero  as  is  immediately  apparent  from  page  256  of 
[Heiskanen  and  Moritz,  1967].  In  particular,  the  elements  of  P'  can  be  expres¬ 
sed  via  a  S.H.  expansion  in  degrees  n+1  through  n'  ,  and  the  elements  of  F"  can 
similarly  be  expressed  in  degrees  n'+l  through  °°.  Thus  any  element  of  P'F"T 
contains  a  sum  of  products  in  which  only  different  degrees  are  present.  Due 
to  the  average  product  (under  the  operator  M)  of  two  Laplace  harmonics  of 
different  degrees  being  zero,  it  follows  that 

M(P‘ F"t}  =  0  ,  (6.33b) 

and,  similarly, 

M{F'F,,T}  =  0  .  (6.33c) 

With  (6.33a,b)  equation  (6.32)  becomes 

M{P'FbT}  =  M{P'F'T}  .  (6.34) 

In  considering  (6.27)  and  (6.28),  we  develop 

M{FbFbT)  =  M{F'F,T}  +  M{F'e' T)  +  M{e'F,T)  +  M{e'e'T}  ,  (6.35) 

where,  similar  to  (6.31)  and  (6.33a), 

M{F'F't}  =  M{F,F,T}  . 

Furthermore, 

M{F'e'T}  =  M{F'F"t}  +  M(F'eT}  i  M{F,F,,T}  +  M(F'  }E(eT}  =  0  , 
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where  use  was  made  of  (6.33c).  Finally,  we  have 

M{e'e'T}  =  M{F"F,,t}  +  M{F"eT}  +  M{eF"T}  +  M{eeT}  s  M{F"F"T}  +  E{eeT)  . 

Upon  collecting  the  last  three  results,  (6.35)  yields 

M{ FbFbT)  =  M{F'F,T}  +  M{F"F,,T}  +  E{eeT}  .  (6.36a) 

This  can  also  be  written  as 

M{FbFbT}  =  M{FFt)  +  E{eeT}  ,  (6.36b) 

due  to  the  fact  that 

M{FFt}  =  M{F'F,t}  +  M{F'F"t}  +  M{F"F,t}  +  M{F" F"T)  =  M{F'F,t}  +  M{F"F"t}  , 

where  (6.33c)  has  been  utilized.  It  is  thus  confirmed  that  the  matrix  in  (6.36a) 
is  the  same  as  that  expressed  by  (6.8). 

In  applying  (6.31),  (6.34)  and  (6.36b)  to  (6.30),  we  have 

M{eVT}  =  M{P'P,t}  -  MfP’F'^A'  -  A,tM{F’P,t}  +  A,t(M{FFt}  +  E{eeT})A'  . 

This  is  written  as 

M{gVt}  =  H  -  M,tA'  -  A,tM'  +  A,t(H  +  z)A*  ,  (6.37) 

where  the  new  notations  are 

M't  =  M{P'F,t}  ,  Hp,E  M{P'P't}  . 

Again  in  analogy  to  Appendix  2,  the  matrix  in  (6.37)  should  have  a  minimum  trace. 
In  particular,  if  one  applies  the  operator  "Tr"  and  differentiates  with  respect 
to  A,T,  the  result  has  to  be  a  zero  matrix: 


3Tr[Hp,]/3A'T  -  23Tr[A,TM']/3A'T  +  3Tr[A'T(H  +  E)A']/3A,T  =  0. 

The  rules  for  differentiation  of  traces  yield 
0  -  2M't  +  2A,t(H  +  E)  «  0  , 
and  thus 

A,t  =  M,t(H  +  z)"1  . 

It  now  follows  from  (6.29)  and  (6.37)  that 

P'  =  M,TH-1Fb  ,  (6.38) 

Cp,  =  Hpl  -  ,  (6.39) 

where  H  has  the  same  meaning  as  in  (6.24).  We  may  add  that  with  the  self- 
evident  changes  in  notations,  (6.36a)  allows  H  to  be  also  expressed  as 

H  =  H'  +  (H"  +  E)  , 


which  corresponds  to  the  splitting  in  (6.26),  in  analogy  to  the  correspondence 
between  (6.23)  and  (6.24).  In  comparing  (6.21)  and  (6.22)  with  (6.38)  and  (6.39), 
we  see  that  the  only  difference  between  the  formulas  pertaining  to  the  actual 
geoid  and  the  smoothed  geoid  resides  in  replacing  MT  and  Hp  in  the  former  by 
M,t  and  ;  this  is  accomplished  simply  by  replacing  "<*>"  (in  practice  1,000, 
etc.)  by  n'  in  conjunction  with  D(^).  As  has  been  already  indicated,  the  matrix 
H  given  in  (6.24)  is  unchanged,  i.e.,  the  computation  proceeds  with 

Besides  the  computational  shortcut,  the  above  procedure  is  compatible  with 
expressing  the  geoid  through  the  S.H.  expansion  to  the  degree  and  order  (n',n'). 
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In  particular,  one  can  utilize  the  predicted  values  P‘  distributed  in  an 
appropriate  (equilateral)  grid  and  compute  the  S.H.  potential  coefficients 

A 

via  integral  formulas.  If  "unsmoothed"  values  P  were  utilized  instead,  the 
smoothing  would  occur  at  the  level  of  these  coefficients.  As  a  final  remark, 
we  observe  that  the  above  derivation  is  general,  in  the  sense  that  the  "un¬ 
smoothed"  formulas  (6.21)  and  (6.22)  follow  from  (6.38)  and  (6.39)  simply  by 
replacing  n'  by  ~  as  a  special  case. 

One  of  the  final  products  we  wish  to  form,  using  the  geoidal  predictions 
obtained  via  the  modified  L.S.  collocation  with  noise  as  has  just  been  described, 
is  represented  by  contour  maps  of  (smoothed)  geoid  undulations  and  gravity 
anomalies  over  the  oceanic  regions  covered  by  satellite  altimetry.  The 
resolution  implied  by  such  maps  corresponds  to  the  density  of  the  prediction 
points  and  to  the  smoothness  of  the  predictions  implied  by  the  degree  n'. 

Clearly,  these  two  characteristics  should  be  consistent  with  each  other.  Thus, 
if  n'  corresponds  to  a  2°  geoidal  resolution,  the  prediction  points  should  form 
a  2°x2°  equilateral  grid,  not  a  coarser  grid.  A  finer  grid  would  not  detract 
from  the  resolution  displayed  by  a  contour  map  but  would  be  uneconomical, 
depending  on  its  interval  and  other  characteristics  (e.g.,  a  l°xi°  geographical 
grid  would  be  wasteful  compared  to  a  l°xi°  equilateral  grid  which,  in  turn, 
would  be  wasteful  compared  to  a  2°x2°  equilateral  grid).  The  appropriate  size 
of  the  grid  can  also  serve  as  a  guidance  when  determining  a  useful  size  of  the 
neighborhood  (around  a  given  prediction  point)  supplying  the  data  values  Fb  in 
(6.38);  for  example,  observation  points  removed  several  grid  intervals  from 
this  prediction  point  would  have  little  effect  on  the  predicted  value  and  could 
be  left  out  of  the  computation. 
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We  now  focus  our  attention  on  developing  the  procedure  for  drawing 
geoidal  and  gravity  anomaly  maps  consistent  with  the  resolution  implied  by  n', 
such  as  a  2°  resolution.  Given  a  2°x2°  equilateral  grid  of  predicted  geoid 
undulations  referring  to  the  (n,n)  field,  various  interpolation  techniques 
could  be  used  to  fill  in  the  undulations  at  intermediate  points  as  needed 
by  a  particular  contour  routine.  For  example,  a  polynomial  interpolation 
could  be  used  for  this  purpose  (some  contour  routines  themselves  use  this  type 
of  interpolation).  The  intermediate,  or  densifying,  points  could  be  chosen 
to  form  a  1°*10  geographical  grid  from  which  a  geoidal  map  could  be  drawn  in 
the  Mercator  projection,  etc.  It  is  understood  that  all  the  predicted  geoid 
undulations  should  be  added  algebraically  to  the  undulations  computed  through 
the  (n,n)  S.H.  expansion  in  order  to  arrive  at  the  final  geoid  undulations 
corresponding  to  the  (n',n')  field. 

A  suitable  interpolation  technique  yielding  the  desired  densifying  grid 
is  the  best  linear  prediction  method,  or  errorless  collocation,  described  in 
Appendix  2.  Indeed,  all  we  wish  to  accomplish  at  this  point  is  to  density  the 
original  (equilateral)  grid  of  predicted  values  considered  as  perfect  in  view  of 
contour  lines.  Quite  naturally,  the  smoothness  in  the  fill ed-in  values  should 
correspond  to  the  smoothness  in  the  original  predictions.  This  requirement 
is  satisfied  if  the  n'  adopted  for  D(ip)  in  the  errorless  collocation  is  the 
same  as  the  n'  used  previously  in  the  modified  L.S.  collocation  with  noise 
when  formulating  the  matrices  M'T  and  Hpl.  Such  a  simple  relationship  between 
the  densifying  method  and  the  original  prediction  method  —  fulfilling  the  need 
for  uniform  resolution  characteristics  --  demonstrates  the  suitability  of  the 
errorless  collocation  for  the  task  at  hand. 
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In  accordance  with  the  above  discussion,  the  following  prediction 
algorithm  can  be  used  for  grid  densifications: 

P‘  =  M*th;;P‘  ,  (6.40) 

C|.  =  H~ ,  -  M-VM*  ,  (6.41) 

where  the  symbol  indicates  the  densified  geoid  undulations.  In  this  case  P' 
in  (6.40)  designates  the  vector  of  predicted  values  obtained  in  (6.38)  and  Hp,  in 
(6.40)  and  (6.41)  designates  the  matrix  which  appeared  in  the  same  notation  on  the 
right-hand  side  of  (6.39).  The  matrices  M,T  and  Hp,  are  formed  according  to  the 
locations  of  the  densifying  points  via  the  covariance  function  D( i^) .  If  a  densi- 
fying  point  coincides  with  an  original  prediction  point,  the  original  predicted 
quantity  is  exactly  reproduced  by  this  method  (see  Appendix  2) ,  which  is  a  desired 
interpolation  characteristic  from  the  standpoint  of  contour  lines. 

But  the  equations  of  the  type  (6.40)  and  (6.41)  can  be  used  for  grid  densi¬ 
fications  in  terms  of  other  geophysical  quantities  as  well.  Assume  first  that 

(6.38)  and  (6.39)  have  been  used  for  smoothed  out  predictions  of  gravity  anomalies, 
in  which  case  the  covariance  function  D(if>)  serving  in  the  formation  of  M'^  and  Hp, 
has  been  replaced  by  the  cross-covariance  function  H(^)  and  the  covariance  function 
C(^),  respectively.  Next  assume,  for  a  moment,  that  the  notations  in  (6.38)  and 

(6.39)  are  retained  also  in  this  new  situation.  It  then  follows  that  the  last 
three  sentences  of  the  preceding  paragraph  can  be  repeated  as  they  stand  also  with 
regard  to  the  densifications  of  gravity  anomalies,  except  that  C (ip)  replaces  D(4>). 

The  values  which  could  serve  in  evaluating  the  effectiveness  and  the 
resolution  power  of  the  final  product  are  the  residuals  computed  with  respect 
to  the  geoidal  contour  map.  One  could  locate  a  specific  observation  point  on 


the  map,  estimate  its  geoid  undulation  referring  to  the  (n,n)  field,  and 
compute  the  difference  between  this  value  and  its  counterpart  in  Fb.  However, 
such  an  estimated  undulation  can  be  found  much  more  easily  upon  using  the  same 
procedure  that  has  led  to  the  construction  of  the  geoidal  contour  lines  them¬ 
selves.  In  particular,  the  formula  (6.40)  can  be  used  in  this  task,  where 
the  observation  points  are  treated  exactly  in  the  same  fashion  as  the  densify- 
ing  points.  In  practice,  only  a  limited  area  associated  with  the  values  P' 
is  considered  around  the  observation  point,  similar  to  the  procedure  mentioned 
earlier  with  regard  to  the  neighborhood  containing  the  observations  Fb.  Quite 
logically,  the  size  of  the  area  used  in  the  residual  computation  should  be 
compatible  with  its  counterpart  in  the  densification  (6.40);  the  approximations 
committed  in  computing  the  residuals  will  thus  be  compatible  with  the  approxima¬ 
tions  embodied  in  the  construction  of  the  contour  lines. 
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6.4  Practical  Considerations 


At  the  outset  we  recall  two  phases  of  the  L.S.  adjustment  used  at  AFGL 
in  satellite  altimetry  reductions:  the  first,  spherical -harmonic  adjustment 
yields  a  smoothed  out  global  geoid,  and  the  second,  point-mass  adjustment 
results  in  a  more  detailed  geoid  on  a  local  or  a  large-scale  basis.  The 
residuals  computed  after  both  adjustment  phases  contain,  in  addition  to  a  ran¬ 
dom  component  (due  to  observational  noise  and  unmodeled  tidal  and  other  effects), 
also  unaccounted  for  short-wavelength  geoidal  information. 

These  residuals  can  be  exploited  in  their  own  right.  For  example,  upon 
using  some  approximate  procedure  they  could  be  translocated  to  form  a  postulated 
grid  and  then  treated  by  various  computationally  efficient  techniques  to  yield 
a  "third -phase"  (detailed)  resolution.  However,  the  collocation  approach 
presented  in  Section  6.3  can  also  lead  to  such  a  third -phase  resolution,  and 
can  be  used  on  a  local,  regional,  or  global  basis.  Furthermore,  it  yields 
error  estimates  associated  with  the  predicted  quantities.  Clearly,  this 
technique  can  be  based  directly  on  the  residuals  from  the  global  adjustment  as 
in  the  second-phase  resolution  which,  so  far,  has  been  accomplished  through 
the  point-mass  adjustment.  The  principle  of  the  collocation  application  remains 
the  same  for  either  choice,  i.e.,  for  the  second-phase  or  the  third-phase 
resolution. 

From  the  conceptual  point  of  view,  the  second-  or  third-phase  collocation 
process  could  be  called  an  "observation  domain"  solution  since  the  output 
quantities  are  of  the  same  kind  as  the  input  quantities, namely  the  geoid  undula¬ 
tions.  With  the  aid  of  the  appropriate  cross-covariance  function,  the  geoid 
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undulations  can  be  subsequently  transformed  into  gravity  anomalies,  etc. 

We  recall  that  *his  approach  does  not  use  intermediate  parameters  such  as  the 
point-mass  magnitudes  computed  as  a  part  of  the  previous  second-phase  L.S. 
adjustment  process. 

Using  the  altimeter  data  entails  some  practical  difficulties  which  could 
hinder  the  efficiency  of  the  collocation  technique  from  the  standpoint  of  both 
the  computer  storage  and  run-time  requirements.  Most  notable  in  this  respect 
is  the  presence  of  a  great  many  observation  points  within  small  oceanic  areas. 

In  addition  to  the  computer  burden,  this  phenomenon  could  also  cause  numerical 
problems  during  the  inversion  of  the  matrix  designed  as  H  in  Section  6.3.  In 
particular,  if  any  two  observation  points  were  located  in  close  proximity  to 
each  other,  H  could  become  ill-conditioned,  especially  if  the  noise  level 
represented  by  the  matrix  E  is  low. 

The  above  difficulties  can  be  circumvented  by  using  the  algorithm  developed 
in  Appendix  3,  which  could  be  viewed  as  a  specific  kind  of  translocation.  It 
performs  a  weighted  averaging  of  the  location  of  two  observation  points  situated 
closer  together  than  a  pre-determined  limit,  as  well  as  the  same  weighted 
averaging  of  the  corresponding  geoidal  residuals  (the  latter  form  Fb  of  Section 
6.3,  i.e.,  the  input  for  the  collocation  predictions).  If  the  resulting 
observation  point  is  too  close  to  some  other,  a  similar  process  again  takes 
place.  However,  instead  of  averaging  the  geographic  coordinates  <p  and  a,  the 
algorithm  gives  directly  the  desired  trigonometric  functions.  Thus,  instead  of 
nine  columns  per  observation  point  only  seven  columns  are  needed  in  a  "storage" 
matrix  whose  rows  correspond  to  the  number  of  the  prediction  points.  The 
locations  of  the  prediction  points  are  usually  stipulated  beforehand;  as  we 
have  seen  earlier,  it  is  advisable  that  they  form  an  equilateral  grid. 
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The  elimination  of  the  above  two  columns  —  which  are  not  needed  in  the 
actual  prediction  process  --  allows  for  the  inclusion  of  more  prediction  points 
and/or  more  observation  points  per  prediction  point  than  would  otherwise  be 
possible  in  a  computer  run.  The  savings  amount  to  about  29%  of  the  above 
"storage"  matrix  which  has  by  far  the  greatest  storage  requirements  of  the 
entire  collocation  application.  The  basic  seven  columns  mentioned  in  the  last 
paragraph  represent  the  following  information:  costj^,  sin^,  cosA^  sinAi, 
cos^ji,  sin^ji  and  ,  where  the  index  J  pertains  to  a  specific  prediction 
point  and  the  index  i  pertains  to  the  observation  point  whose  input  geoid 
undulation  (more  precisely,  the  part  unmodeled  by  the  previous  adjustment  plus 
the  observational  noise)  is  N.. 

The  observation  points  i  are  arranged  according  to  their  angular  distance 
(<J» )  from  the  prediction  point  J,  so  that  the  point  No.  1  is  closer  to  J  than  is 
the  point  No.  2,  etc.  If  neither  of  the  considered  points  has  been  averaged 
before,  the  respective  weight  factors  are  taken  as  unity.  If  a  point  has  been 
averaged  once  before  and  if  another  averaging  is  to  take  place,  its  weight 
factor  is  2.  Except  for  the  weight  factor,  the  weight  is  in  reverse  proportion 
to  the  (angular)  distance  between  the  observation  point  and  its  prediction 
point.  This  importance  given  to  the  distance  is  rooted  in  the  consideration  of 
the  prediction  formula  where  the  contribution  of  individual  observation  points 
diminishes  quite  rapidly  with  distance  from  the  prediction  point,  reflecting 
the  behavior  of  the  covariance  function. 

As  has  been  already  pointed  out  in  Section  6.3  the  basic,  or  original, 
predictions  are  considered  to  form  an  equilateral  grid.  Such  a  distribution 
gives  rise  to  a  uniform  resolution  depending  on  the  grid  interval  (the  degree 
n'  defined  in  the  same  section  is  assumed  to  be  compatible  with  this  interval). 


Accordingly,  a  unit  sphere  representing  the  globe  is  thought  of  as  subdivided 
into  compartments,  each  characterized  by  the  appropriate  angular  intervals 
A<}>  and  AX,  and  each  associated  with  one  centrally  located  prediction  point. 

The  value  of  A<j>  is  independent  of  position,  but  the  value  of  AX  varies  as 
A<t>/ co sd> ,  where  <p  is  the  geocentric  latitude.  The  poles  would  be  an  exception 
to  this  rule  (the  compartments  would  become  spherical  caps),  but  they  are 
omitted  from  this  development.  In  the  case  of  SEASAT,  for  example,  the  ground 
tracks  reach  only  the  latitudes  +.72 0  beyond  which  no  predictions  are  made. 

The  continental  areas  with  no  altimeter  measurements  are  likewise  ignored  in 
the  current  application. 

Although  in  theory  all  the  observation  points  should  be  involved  in  any 
prediction,  in  practice  such  a  procedure  would  be  extremely  wasteful.  Due  to 
the  nature  of  the  geoidal  covariance  function,  the  observation  points  represented 
by  the  vector  Fb  can  be  restricted  to  a  relatively  small  neighborhood  of  the 
pertinent  prediction  point  (called  J).  Clearly,  all  the  observation  points 
from  the  compartment  represented  by  J  should  be  involved.  However,  the  question 
of  how  far  beyond  this  compartment  the  search  for  observation  points  should  be 
extended  is  a  matter  of  judgment  and  experience,  depending  to  a  large  extent 
on  the  cmputer  characteristics  and  limitations. 

Practical  considerations  coupled  with  numerical  evaluations  have  indicated 
that  the  inclusion  of  observation  points  located  within  the  spherical  cap  that 
just  covers  the  corners  of  the  compartment  J  is  reasonable.  This  provides  for 
a  modest  measure  of  continuity  between  neighboring  predictions  because  of  the 
overlap  in  peripheral  observations.  The  radius  of  such  a  cap,  centered  at  J, 
is  approximately  0.75a$.  Even  so,  the  number  of  observation  points  involved 
in  one  prediction  may  sometimes  be  prohibitive.  This  occurs  in  spite  of  the 


input  data  being  limited  to  0.5°-intervals  along  satellite  arcs,  and  is  due 
to  a  substantial  increase  in  the  number  of  observations  brought  about  by  the 
presence  of  repeating  tracks.  It  then  becomes  necessary  to  stipulate  a 
maximum  allowable  number  of  observation  points  and  discard  the  points  farthest 
away  from  J  until  this  limit  is  satisfied.  For  this  reason  it  is  useful  to 
have  the  observation  points  ordered  according  to  their  angular  distances  (^) 
from  J. 

The  number  of  observation  points  within  the  spherical  cap  introduced 
above  depends,  of  course,  on  the  desired  resolution.  In  a  2°  resolution  (the 
grid  size  corresponds  to  A4>=2 0 )  the  number  of  observations  often  reaches 
beyond  50.  In  a  1°  resolution  this  number  is  likely  to  decrease  four-fold. 

If  the  limit  in  the  former  case  is  set  to  50-60,  most  of  the  observations 
within  the  pertinent  cap  will  be  accepted.  On  the  other  hand,  the  limit  of 
20-30  in  the  latter  case  will  often  result  in  all  of  the  observations  being 
accepted.  These  limits  serve  merely  as  practical  guidelines.  Clearly,  the 
final  number  of  observation  points  is  also  reduced  by  the  special  translocation 
(averaging)  technique  described  in  the  first  part  of  this  section.  In  this 
respect,  the  minimum  separation  limit  of  0.2A<f>  has  proven  useful. 

With  regard  to  the  densifying  points,  it  is  suggested  that  the  spherical 
cap  centered  at  any  such  point  be  enlarged  to  over  1A$  (e.g.,  1 . 125Ad>  has  been 
used  satisfactorily).  This  suggestion  stems  from  the  fact  that  only  relatively 
few  prediction  points,  represented  by  the  vector  P1  in  (6.40),  are  located 
around  any  densifying  point.  Similar  to  an  earlier  statement,  the  points 
where  the  "residuals"  are  sought  should  again  be  treated  in  the  manner  of 
densifying  points. 
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Orbital  considerations  dictate  that  the  finest  grid  of  prediction  points 
to  be  used  in  conjunction  with  SEASAT  observations  is  a  l°xl°  equilateral  grid. 
This  limitation  is  imputable  to  the  ground  tracks  intersecting  in  a  very 
approximate  equilateral  grid  with  about  1°  on  the  side  over  much  of  the  oceanic 
surface.  The  contour  maps  could  thus  represent  a  1°  resolution  at  best,  which 
is  understood  as  expressing  the  features  down  to  and  including  those  whose 
shortest  half-wavelength  is  1°.  Since  such  a  resolution  corresponds  approximate¬ 
ly  to  a  (180,180)  set  of  S.H.  potential  coefficients,  the  truncation  n'  used  in 
Section  6.3  in  conjunction  with  the  covariance  and  cross-covariance  functions 
is  approximately  n'=180. 

Finally,  we  note  that  the  predictions  described  herein  can  be  considered 
as  "best",  but  only  when  viewed  as  a  whole,  i.e.,  on  the  global  scale.  This 
stems  from  the  global  character  of  the  covariance  and  cross-covariance  functions 
characteristic  of  the  M-averaging  as  presented  in  Appendix  2  and  as  used  in 
Section  6.3.  If  a  similar  procedure  should  be  applied  on  a  regional  basis, 
regional  covariance  and  cross-covariance  functions  as  well  as  regional  averaging 
would  have  to  be  employed  in  order  to  yield  some  "best"  predictions  on  a 
smaller  scale.  If  such  regional  functions  —  or  any  other  functions  different 
from  those  in  Appendix  2  —  were  used  for  predictions  covering  the  entire  globe, 
the  quality  of  these  predictions  as  a  whole  would  necessarily  be  worse  than  the 
quality  of  the  predictions  dealt  with  herein. 


-79- 


7.  CONCLUSIONS 


The  topics  addressed  in  this  report  have  separate  characteristics.  They 
are  treated  in  Chapters  2-6  which  are  essentially  self-contained.  Accordingly, 
only  a  brief  summary  of  most  topics  is  presented  below. 

In  Chapter  2,  a  classical  oceanographic  procedure  is  justified  by  showing 
that  except  for  a  factor,  the  pressure  gradient  along  equi potential  surfaces 
is  interchangeable  in  practice  with  the  geopotential  gradient  along  isobaric 
surfaces.  The  use  of  the  normal  coordinate  system  allows  the  assessment  of 
the  approximations  involved  in  such  a  practical  formulation,  as  well  as  the 
development  of  the  rigorous  formulation  with  an  additional  factor  containing 
dW/g,  where  W  is  the  geopotential  and  g  is  the  gravity.  Except  for  the  geo¬ 
potential  considered  as  increasing  downward  according  to  the  geodetic  conven¬ 
tion,  the  tensor  notations  and  concepts  of  [Hotine,  1969]  are  adhered  to. 

Chapter  3  takes  advantage  of  the  formula  for  the  bottom  tide  as  developed 
in  [Schwiderski ,  1980]  and  in  other  reports  by  the  same  author.  Although  these 
reports  do  not  deal  with  altimeter  adjustment  per  se,  an  improvement  of  the 
altimeter  model  with  tidal  effects  is  achieved  through  the  replacement  of  the 
approximate  value  of  the  ocean  bottom  deformation  due  to  ocean  tidal  loading 
by  a  more  rigorous  formulation.  As  an  additional  advantage,  this  formulation 
can  be  used  in  conjunction  with  all  of  the  tidal  constituents. 

In  Chapter  4,  an  economical  yet  quite  rigorous  weighting  scheme  is 
developed  for  altimeter  observations  in  the  first-phase  adjustment.  It  is 
based  on  the  first-order  autoregressive  process  of  [Brown  and  Trotter,  1969] 
applied  to  the  values  obtained  with  the  geoidal  covariance  function.  Its 
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main  outcome  is  represented  by  a  tri-diagonal  weight  matrix  which,  when 
compared  to  the  rigorous  case  of  a  fully  populated  weight  matrix,  offers 
almost  the  same  level  of  computer  economies  as  a  diagonal  weight  matrix. 
However,  the  errors  associated  with  a  diagonal  matrix  are  greatly  reduced 
in  this  approach,  perhaps  by  60%. 

Chapter  5  extends  the  concept  of  second-phase  adjustment  of  satellite 
altimetry  by  incorporating  chosen  tidal  parameters  in  the  point-mass  model. 
Although  the  point  masses  in  a  given  ocean  basin  receive  in  general  different 
adjustment  corrections,  it  is  not  so  for  the  tidal  parameters  considered 
common  to  the  whole  area.  These  parameters  are  now  associated  with  the 
border  of  the  matrix  of  normal  equations.  Accordingly,  whereas  the  "modified 
Choleski  algorithm"  of  [Blaha,  1983  ]  containing  the  point-mass  parameters 
alone  corresponds  to  a  banded  system  of  normal  equations,  the  new  "extended 
modified  Choleski  algorithm"  results  in  a  banded-bordered  system  of  normal 
equations.  In  the  closing  part  of  Chapter  5  it  is  confirmed  that  the  former 
is  a  special  case  of  the  latter.  In  either  case,  the  presented  algorithm  is 
computationally  efficient. 

Finally,  Chapter  6  offers  an  alternative  to  a  second-  or  third-phase 
simultaneous  least-squares  adjustment  process.  In  particular.  Section  6.3 
contains  an  approach  based  on  the  development  of  the  collocation  theory  as 
presented  in[Heiskanen  and  Moritz,  1967]  and  [Moritz,  1980].  Both  the 
signal  and  the  noise  are  considered  to  be  present  at  any  observation  point, 
and  are  associated  with  the  averaging  operators  M  and  E,  respectively.  The 
derivations  proceed  strictly  in  matrix  notations.  It  is  shown  that  if  the 
predicted  values  of  geophysical  quantities  should  he  desired  in  a  smoothed 


out  version,  the  standard  collocation  approach  is  unchanged  except  that 
during  the  formation  of  the  first  matrix  in  the  prediction  formula,  the 
covariance  (or  cross-covariance)  function  is  no  longer  utilized  in  conjunc¬ 
tion  with  the  expansion  degree  but  only  in  conjunction  with  the  degree 
n'  corresponding  to  the  required  smoothness.  As  n'  increases,  the  standard 
"unsmoothed"  collocation  version  is  seen  to  follow  from  its  smoothed  out 
counterpart  as  a  special  case. 

A  somewhat  more  detailed  sunmary  of  the  most  important  outcome  of 
Section  6.3  is  deemed  useful  (without  the  re-introduction  of  the  notations 
utilized  therein).  With 


H  =  H  +  Z  , 


equations  (6.21)  and  (6.22)  are  recapitulated  as 


P  =  MWb  , 


Cp  =  Hp  -  MTH_1M  . 


The  matrices  M  ,  H  and  Hp  are  formed  using  the  covariance  function  D(^)  cor¬ 
responding  to  (the  signal  part  of)  Fb.  The  latter  consists  here  of  geoid 
undulations,  but  this  case  can  be  generalized  with  minimal  and  self-evident 
changes. 

With  regard  to  geophysical  quantities  other  than  those  represented  by 
Fb,  one  can  go  over  the  derivations  leading  to  the  above  formulas  and  add 
the  symbol  to  all  the  vectors  and  matrices  except  Fb,  H,  Z  (and  W  if  it 
should  be  used).  The  result  is 
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P  =  , 

C~  =  H~  -  . 

P  P 

In  the  current  application  the  desired  quantities  are  gravity  anomalies;  the 
matrix  MT  is  thus  formed  using  the  cross- covariance  function  H(i/>)  and  H~  is 
formed  using  the  covariance  function 

Equations  (6.38)  and  (6.39)  pertaining  to  the  smoothed  out  predictions 
are  recapitulated  as 

P'  =  M,TH_1Fb  , 

C$.  =  Hp<  -  M' TH-1M'  . 

The  matrices  M,T  and  Hp,  are  computed  as  their  unprimed  counterparts  except 
that  the  degree  "<=°"  in  the  formation  of  D( ip)  is  replaced  by  n'.  As  n'  in¬ 
creases  to  infinity,  these  formulas  become  their  unprimed  counterparts. 

If  one  seeks  the  smoothed  out  predictions  of  geophysical  quantities  other 
than  those  represented  by  Fb,  in  analogy  to  the  case  already  treated  one  can 
go  over  the  derivations  and  add  the  symbol  to  all  the  vectors  and  matrices 
except  Fb,  H,  E  (and  W  if  it  should  be  used).  The  result  is 

P'  =>  M'TH-1Fb  , 

C~,  =  H~,  -  . 

The  matrices  M'T  and  H~  are  computed  as  their  unprimed  counterparts  except 
that  the  degree  is  again  replaced  by  n\  As  n'  increases  to  infinity. 
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these  formulas  become  their  unprimed  counterparts.  One  notices  that  the 
vector  H_1FbEW  is  unchanged  in  all  four  prediction  formulas  above. 

Section  6.3  also  demonstrates  that  the  densifications  of  the  predicted 
geophysical  quantities,  needed  especially  in  view  of  contour  maps,  can 
proceed  with  advantage  upon  using  the  errorless  collocation.  The  effective¬ 
ness  of  the  presented  approach  in  its  entirety  can  be  evaluated  with  the  aid 
of  the  residuals  (at  observation  points)  computed  with  respect  to  the  geoidal 
contour  map.  Here  the  errorless  collocation  represents  not  only  a  useful  but 
also  a  consistent  tool,  due  to  the  fact  that  it  has  lead  to  the  construction 
of  the  geoidal  contour  map  itself. 
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APPENDIX  1 


COMPUTER  PROGRAM  FOR  GEOIDAL  VARIANCE-COVARIANCES 
DUE  TO  TRUNCATIONS 

IN  THE  UNDERLYING  SPHERICAL-HARMONIC  MODEL 


PROGRAM  CNTRL  (INPUT,  OUTPUT,  TAPE5= INPUT,  TAPE6=OUTPUT) 

COMPUTE  VARIANCE-COVARIANCES  (IN  METERS  SQUARE)  DUE  TO  TRUNCATIONS 
OF  THE  SPHERICAL-HARMONIC  MODEL  FOR  GEO ID  UNDULATIONS 

COMMON/ORDER/ IORDER 
COMMON/ COSPSI/2ED 
C0MM0N/PLEG/PZ(1000) 

DIMENSION  DEGVAR(IOOO) 

DIMENSION  PS I( 30) 

DIMENSION  N (20) 

DIMENSION  SM(20,30) 

DATA  ZER,  RAD/O. ,0.01745329252/ 

LOGIN  =  5 
LGOUT  =  6 

WRITE(LG0UT,5) 

IZ  =  0 

READ(L0GIN,15)  NPROB 

1  IZ  =  IZ  +  1 
WRITE (LG0UT,25)  IZ 
READ(L0GIN,15)  IORDER 

IF( IORDER. GT. 1000)  IORDER  =  1000 
..RITE(LG0UT,35)  IORDER 
READ(LOGIN ,15)  NPSI 
WRITE (LGOUT, 45)  NPSI 
READ(L0GIN,55) (PSI ( I) ,1  =  1 .NPSI ) 

WRITE(LG0UT,55)(PSI(I),I  =  l.NPSI) 

DO  2  I  =  l.NPSI 

2  PSI(I)  =  RAD*PSI( I) 

READ(L0GIN,15)  NN 
WRITE (LGOUT, 65)  NN 
READ(L0GIN,15)(N(J),U  =  1,NN) 

WRITE(LG0UT,15)(N(J),J  =  1,NN) 

DEGVAR(l)  =  ZER 
DEGVAR(2)  =  ZER 
DO  3  K=3, IORDER 
XKM1  =  K  -  1 
XKM2  -  K  -  2 
XKP24  =  K  +  24 
KP2  =  K  +  2 

DEGVAR(K)  =  (0.9996 17**KP2 )*17981 . / ( XKM1*XKM2*XKP24 ) 

3  CONTINUE 
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ZED  =  COS(PSKD) 

CALL  LEGEN 
K  =  IORDER  +  1 
SUM  =  ZER 
IPART  =  NN  +  1 

6  CONTINUE 

IPART  =  IPART  -  1 
NPART  =  N( IPART) 

7  CONTINUE 
K  5  K  -  1 

IF(K.EQ. NPART)  GO  TO  8 
SUM  =  SUM  +  DEGVAR(K)*PZ(K) 
GO  TO  7 

8  CONTINUE 

SM ( IPART, I )  =  SUM 
K  =  K  +  1 

IF( IPART. GT.l)  GO  TO  6 
4  CONTINUE 


DO  9  IPART  =  1 ,NN 
WRITE(LGOUT,75)  N( IPART) 
WRITE(LG0UT,55) (SM( IPART , I ) , I 
9  CONTINUE 

IF( IZ.LT.NPROB)  GO  TO  1 
STOP 


=  l.NPSI) 


5  FORMAT(///) 

15  FORMAT( 1615) 

25  F0RMAT(////,15X,11HPR0BLEM  NO., 13) 

35  FORMAT (//,5X,46H INFLUENCE  CF  NEGLECTED  DEGREES  SUMMED  UP  TO  N=,I3) 

45  FORMAT(//,5X,3HTHE,I3,27HVALUES  OF  PSI  (DEGREES)  ARE) 

55  FORMAT (8F 10. 3) 

65  FORMAT ( //  ,5X , 3HTHE  ,13, 38HVALUES  OF  INDIVIDUAL  TRUNCATIONS  N  ARE) 

75  F0RMAT(///,2X,50HEFECT  (FOR  GIVEN  ANGLES  PSI)  OF  THE  TRUNCATION  N  =,I3,/) 
END 


SUBROUTINE  LEGEN 
COMMON/ORDER/ IORDER 
COMMON/ COSPS I/ZED 
C0MM0N/PLEG/PZ(1000) 
DATA  ONE,  TWO/1. ,2./ 
PZ(1)  =  ZED 
PM2  =  ONE 
PM1  =  ZED 

DO  100  N  =  2, IORDER 
XN  =  N 

VAR1  =  ONE/XN 
VAR2  =  TWO*XN-ONE 


B  =  VAR1*(ZED*VAR2*PM1  -  VAR3*PM2) 
PZ(N)  =  B 
PM2  =  PM1 
PM1  =  B 
100  CONTINUE 
RETURN 
END 


\ 

\ 
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APPENDIX  2 

BEST  LINEAR  PREDICTION  METHOD  ON  A  GLOBAL  SCALE 


The  development  in  this  analysis  follows  closely  that  of  [Heiskanen 
and  Moritz,  1967],  except  that  the  matrix  conventions  are  used  throughout. 

It  is  based  on  the  notion  of  covariance  function  as  presented  in  Chapter  7 
of  this  reference,  which  plays  a  fundamental  role  in  the  process  of  predic¬ 
tions  (interpolations  or  extrapolations).  Two  remarks  with  regard  to  the 
adopted  methodology  are  in  order.  First,  the  symbol  E  for  mathematical 
expectation  is  replaced  by  M,  indicating  averages  over  the  unit  sphere  as 
explained  in  the  above  reference.  And  second,  the  errors  of  the  measurements 
themselves  are  considered  to  be  zero.  The  errors,  averages,  etc.,  are 
associated  exclusively  with  errors  of  prediction.  Therefore,  the  "observa¬ 
tions"  play  the  role  of  errorless  theoretical  quantities,  which  means 
that  they  equal  their  expected  Vc  jes.  The  predictions  computed  at  points 
coinciding  with  "observation  points"  will  be  seen  to  be  equal  to  the  original 
observations  and  their  error  estimates  (defined  in  terms  of  M-averages  and  not 
E-averages)  will  have  values  of  exactly  zero. 

One  can  express  all  the  observations,  whose  number  is  denoted  by  n,  as 
the  (column)  vector  F  with  n  elements.  For  the  present  purpose,  they  will  be 
considered  as  consisting  of  geoid  undulations  (little  change  in  the  derivations 
would  occur  if  they  consisted  of  gravity  anomalies  or  other  functions  of  the 
geopotential  instead).  The  quantities  to  be  predicted  will  be  grouped  in  the 

A 

vector  P  which  represents  an  estimate  of  the  vector  P,  containing  unknown 
values  of  some  function  at  arbitrarily  chosen  locations.  This  function  may 


or  may  not  be  the  same  as  the  one  observed,  i.e.,  the  one  whose  selected 
values  are  depicted  by  the  elements  of  F.  In  fact,  two  kinds  of  functions 
will  be  considered  herein  in  connection  with  the  vector  P:  1)  geoid  un¬ 
dulations  as  above,  and  2)  gravity  anomalies.  The  prediction  of  P  will  be 
made  possible  through  the  known  behavior  of  these  geophysical  quantities 
embodied  in  the  covariance  function  (if  P  and  F  are  of  the  same  kind)  or  in 
the  cross-covariance  function  (if  P  and  F  are  of  a  different  kind).  It  is 
to  be  emphasized  that  we  are  concerned  here  with  the  best  linear  prediction, 

A 

in  the  sense  that  P  consists  of  linear  functions  of  the  elements  in  F,  and 
that  the  difference  between  P  and  P  is  minimized  in  a  certain  sense  (upon 
using  the  operator  M). 

The  pertinent  covariance  and  cross-covariance  functions  are  now  described, 
starting  with  the  surface  spherical -harmonic  expansion  of  gravity  anomalies 
(Ag): 

00  n  _  _  _ 

Ag  =  l  y  (c  cos  mA  +  d  sin  mA)P  (sin$)  ,  (A2.1) 

3  L.  Lr,  nm  nm  nm 

n=2  m=0 


valid  for  the  reference  ellipsoid  of  the  same  mass  and  the  same  potential  as  the 
geoid;  see  [Heiskanen  and  Moritz,  1967],  page  252  or  the  equation  (7-13),  etc. 

In  this  formula,  <p  and  A  are  the  geocentric  latitude  and  longitude  of  the 
point  associated  with  Ag.  The  overbars  indicate  that  the  normalization  has 
taken  place  (thus  Pnm  are  the  normalized  associated  Legendre  functions).  The 
relation  between  the  coefficients  in  (A2.1)  and  the  potential  coefficients 
(C's  and  S's  if  normalized,  C's  and  S's  if  conventional)  is  the  following: 


(A2.2a) 
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and  where 


AC  „  ...  the  correction  to  the  reference  C*  in  order  to  obtain  C  , 
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..  all  the  other  coefficients  , 


G  ...  the  average  value  of  gravity. 


It  can  be  shown  that  the  covariance  function  for  gravity  anomalies  is 


(A2.3) 


C(i|>)  s  M{Ag  Ag‘ }  -  l  cn  Pn(cosij;)  , 

n=2 


where  c^  ,  the  degree  variance  for  gravity  anomalies,  is 

c  e  M { Ag 2 >  =  l  ( c2  +  d2  )  =  G2 (n  -  l)2  l  (AC2  +  AS2  )  , 

n  3n  nm  nm  t(.  nm  nm 
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and  t|>  is  the  spherical  distance  between  the  points  associated  with  Ag  and  Ag1. 


Due  to  Pn(cos  0)  e  1  ,  one  also  has 


C(0)  e  M{ Ag2 }  =  l  c  =  l  M{Ag2} 

n=2  n  n=2  n 


(A2.31) 


f 


Since  usually  C2Q=C*0  and  thus  AC2Q=0,  while  ASn0HO  for  any  n,  then 
c2=0  and  the  summations  in  (A2.3)  and  (A2.3‘)  as  well  as  in  similar  expres¬ 
sions  below  can  start  at  n=3.  These  formulas  therefore  pertain  to  the 
"bandwidth"  between  n=3  and  n=°°.  If  the  reference  gravity  field  were  expres¬ 
sed  through  the  degree  (instead  of  2  as  is  the  case  with  the  ellipsoidal 
field)  the  summations  would  start  with  n^l  (instead  of  3).  And  if  the  com¬ 
plete  field  did  not  go  beyond  the  degree  n=n2  the  summations  would  end  with 
n2  (instead  of  °°). 

Under  the  same  assumptions  (equal  mass  and  potential)  as  used  in  (A2.1), 
the  geoid  undulations  (N)  can  be  written  as 

N  =  (R/G)  l  [l/(n  -  1)]  l  (c  cos  mX  +  cf^  sin  mX^Jsint}.)  ,  (A2.4) 

n=2  m=0 

where  R  is  the  earth's  mean  radius.  The  covariance  function  for  geoid  undulations 
can  similarly  be  derived  as 

oo 

D(i)i)  =  M{N  N' }  =  l  dn  Pn(cos<|>)  ,  (A2.5) 

n=2 

where  dn  ,  the  degree  variance  for  geoid  undulations,  is 

d  =  M{N2}  =  (R/[G(n-  l)]}2  c  =  R2  \  (AC2  +  AS2  )  ; 

n  n  n  nm  nm 

from  (A2.5)  it  follows  that 

0(0)  H  M(N2)  =  1  dn  -=  I  M(N^>  .  (A2. 5* ) 

n=2  n=2 
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The  cross-covariance  function  for  gravity  anomalies  and  geoid  undulations 


can  be  obtained  as 

ao 

H(i|>)  =  M{Ag  N'}  e  M{N  Ag ' }  =  7  h  P  (cos-*)  ,  (A2.6) 

n=2  n  n 

where 

h  =  {R/[G(n- l)]}c  =  [G(n  -  1)/R]d  =  RG(n-l)  ?  (AC2  +  AS2  )  , 

n  J  n  n  nm  nm 

m=0 

h  =  M{Ag  N  }  e  M{N  Ag  }  ; 

n  3n  n  n  3n 


one  also  has 


H( 0)  «  [  h  , 

n=2 


(A2.61 ) 


or 

00  CO 

H ( 0 )  =  M{Ag  N}  e  M{N  Ag}  =  l  M{Ag  N  }  E  £  M{N  Ag  }  . 

n=2  n  n  n=2  n  n 

A 

Attention  may  now  be  focused  on  finding  the  best  linear  estimates  P  , 

P  =  AtF  .  (A2.7) 

(This  formulation  resembles  superficially  the  usual  prediction  formula  in 
stochastic  processes;  however,  the  latter  notion  involves  the  expectation 
operator  E  and  permits  ZF  t  0  ,  unlike  the  present  approach.)  The  matrix  AT 
contains  as  yet  unknown  constants.  The  problem  of  finding  AT  is  addressed 
through  the  M  operator.  We  first  form  a  "dispersion  matrix",  where  e=P-P, 
as  follows: 

ecT  e  (P- AtF)(P-AtF)t  =  PPT-PFTA-ATFPT+ATFFTA  ;  (A2.8) 
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its  average  over  the  unit  sphere  is 


M{eeT}  =  M{PPt}  -  M{PFT}A  -  ATM{FPT}  +  ATM{FFT}A  .  (A2.9) 

Although  PPT,  PFT,  FPT  are  unknown,  the  global  means  of  these  values,  M{PPT}, 
M{PFt},  M{FPt}  are  known  from  the  variance  (and  cross-covariance)  function. 

Next  we  stipulate  that  the  matrix  M{eeT}  should  have  a  minimum  trace, 
thus  applying  the  operator  "Tr"  and  differentiating  with  respect  to  AT;  the 
result  has  to  equal  the  zero  matrix,  namely 

9Tr[M{PPT}]/9AT  -  2  9Tr[ATM{FPT}]/9AT  +  9Tr[ATM{FFT}A]/9AT  =  0  , 

where  use  was  made  of  Tr[B] =  Tr[BT].  Applying  the  rules  for  differentiation 
of  traces  (see  e.g.  [Blaha,  1976],  pages  12  and  13)  one  obtains 

0  -  2  M{PFt}  +  2  ATM{FFT}  =  0  , 

and  thus 

At  =  M{PFT}[M{FFT}]_1  .  (A2.10) 

The  estimate  P  is  then  given  from  (A2.7)  as 

P  =  M{PFt}[M{FFt}]“1F  .  (A2.ll) 

If  (A2.10)  is  substituted  into  (A2.9),  the  second,  third  and  fourth  matrices 
on  the  right-hand  side  become  equal  and  the  result  is 

M{eeT}  =  M{PPt)  -  H{PFT}[M{FFT}]-1M{FPT}  .  (A2.12) 
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The  relations  (A2. 10-12)  are  now  written  as 

at  =  mV1  , 
p  =  mth"1f  =  mtw  , 

M{ecT}  »  Hp  -  MtH-1M  , 

with  the  use  of  the  following  notations: 

W  =  H_1F  ; 

H  =  M{FFt}  ; 

Mt  =  M{PFt}  ,  Hp  £  M{PPT)  . 

The  matrices  M,  H  should  not  be  confused  with  the  operator  M{ }  and  the  cross¬ 
covariance  function  respectively.  All  of  the  above  notations  may  be 

used  when  the  functions  in  P  and  in  F  are  of  the  same  kind  (here  geoid  undula¬ 
tions). 

When  these  functions  are  of  a  different  kind  (gravity  anomalies  associated 
with  P,  geoid  undulations  associated  with  F  as  above),  all  the  matrices  and 
vectors  may  be  identified  by  the  symbol  except  F,  H,  W  which  remain  the 
same.  We  can  thus  write 

at  =  mV1  ,  (a  2.10") 

A 

P  =  mV  F  =  MTW  ,  (A2.ll") 

M{ceT}  =  Hg  -  MtH~ .  (A2.12") 


(A2.10*) 
(A2. 11* ) 
(A2.121 ) 
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The  formula  (A2.ll")  is  interesting  in  that  only  the  matrix  MT  changes  (it 
is  now  computed  using  the  cross-covariance  function),  while  the  vector 
W  s  H_1F  is  the  same  whether  P  or  P  are  sought.  This  procedure  can  be  referred 
to  as  a  best  prediction  of  gravity  anomalies  from  geoid  undulations  (errorless 
values  as  stipulated). 

The  matrix  AT  (or  AT)  in  all  previous  formulas  was  found  subject  to  the 
condition  that  the  M-average  of  all  the  diagonal  elements  in  eeT  be  simultane¬ 
ously  as  small  as  possible.  In  considering  only  one  prediction  point,  P,  P,  e 
comprise  each  only  one  element  and  A  ,  M{PF  }  in  (A2.9)  consist  each  of  only 
one  row  as  a  special  case.  The  derivation  would  then  become  somewhat  simpler 
in  that  the  operator  Tr  would  be  omitted.  In  fact,  after  the  transposition  of 
the  third  term  (a  lxl  matrix)  on  the  right-hand  side  of  (A2.9)  and  the  dif¬ 
ferentiation  of  the  whole  expression  with  respect  to  the  vector  A,  one  would 
recover  the  equation  preceding  (A2.10).  The  row  vector  AT  in  (A2.10),  etc., 
would  now  correspond  to  one  minimum  element  M{e2}. 

The  matrices  H,  MT  and  MT  will  now  be  expressed  more  explicitly.  For 
the  sake  of  clarity,  the  observation  points  will  be  identified  with  numbers 
or  small-letter  indices  while  the  prediction  points  will  be  identified  with 
capital -letter  indices.  The  matrix  H  can  be  written  as 


(A2.13) 


H  =  M{FFt}  = 


Dq  D12  ...  Dln 


...a 


2n 


symm. 
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where  D  s  D(0)  and  D.  .  =D(i|). .),  and  where  n,  the  number  of  observation 
points,  need  not  be  confused  with  the  same  symbol  used  earlier  to  denote  the 
degree  in  a  spherical -harmonic  expansion.  The  matrix  MT  is  represented  as 
fol 1 ows : 


Mt  s  M{PFt}  = 


dt1  d  ...  d, 

J1  J2  Jn 


Dt1  Dr0  ...  Dt 

LI  L2  Ln 


(A2.14) 


while  the  matrix  MT  is 


mt  =  M{PFT}  = 


HJ1  HJ2  HJn 


H  H.  ...  H 

LI  L2  Ln 


(A2.15) 


The  matrix  Hp ? M{PPT)  in  (A2.121)  has  a  structure  similar  to  (A2.13)  above, 
except  that  i,j  now  correspond  to  the  predicted  points.  And  the  matrix 
Hp=M{PPT}  in  (A2.12")  would  be  constructed  as  Hp  (i.e.,  using  predicted  points), 
except  that  the  covariance  function  "D"  would  be  replaced  by  the  covariance 
function  "C"  seen  earlier. 

Suppose  now  that  the  point  J  in  (A2.14)  coincides  with  one  point  of  F, 
e.g.  with  the  i-th  point.  The  row  corresponding  to  this  point  in  MT  then 
becomes  [D.,  D,„  ...  0  ...  D.  ],  identical  to  the  i-th  row  in  H.  When  the 

il  i2  o  in 
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multiplication  in  (A2.10)  or  (A2.101)  is  performed,  the  pertinent  row  in  AT 
contains  zeros  everywhere  except  at  its  i-th  place  where  it  contains  +1.  One 
thus  has  for  the  prediction  (A2.ll)  or  (A2.ll'): 


i.e.,  the  prediction  coincides  with  the  (errorless)  observation.  The  same 
argument  may  be  repeated  when  all  the  predicted  points  coincide  with  the  n 
observation  points,  in  which  case 

MT  =  H  , 

P  =  MtH_1F  s  F  , 

so  that  all  the  predictions  coincide  with  the  observations. 

The  situation  with  respect  to  M{eeT)  in  (A2.12)  or  (A2.12')  is  examined 
next.  The  diagonal  of  Hp  contains  Dq  values  everywhere.  When  the  point  J 
coincides  with  the  i-th  point  in  F,  the  pertinent  row  in  MTH-1  again  has  zeros 
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everywhere,  except  for  +1  at  its  i-th  place.  The  corresponding  column  (for 

CT 

D  D  ...  D  ...  D  "I  seen  before  in  the 
il  i2  o  in-* 

transposed  form;  the  element  DQ  appears  at  the  i-th  place.  Thus,  the  element 
of  mV  LM  at  the  location  (i,i)  is  Dq  from  which  it  follows  that 

=  0o  -  0o  =  0  . 

If  again  all  the  predicted  points  coincide  with  the  n  observation  points,  we 
see  immediately  that 

Hp  =  H  , 


and,  as  above, 

T 

M  =  H  ; 

accordingly, 

M(ceT}  =  H  -  HH_1HT  h  0  . 

Thus  the  average  dispersion  matrix  is  zero,  which  could  be  loosely  termed  as 
"perfect  prediction".  It  is  prudent  to  refrain  from  calling  the  diagonal 
elements  in  M{eeT}  as  variances  and  the  off-diagonal  elements  as  covariances, 
since  in  statistics  such  a  terminology  is  customarily  associated  with  random 
variables  (rather  than  with  some  theoretical  function  values)  and  with  the 
operator  E. 

In  summary,  the  theoretical  (errorless)  values  of  some  function  (here 
geoid  undulations)  at  several  discrete  points  are  used  to  make  a  linear 
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prediction  of  the  same  function,  or  a  different  function  (here  gravity 
anomalies),  at  arbitrarily  chosen  locations.  This  is  possible  solely  due 
to  theknowledge  of  behavior  of  these  geophysical  quantities,  embodied  in 
the  covariance  and  cross-covariance  functions;  these  functions  are  obtained 
through  a  stationary  process  on  the  (unit)  sphere  and  are  ultimately  given 
in  terms  of  the  spherical -harmonic  (S.H.)  potential  coefficients.  An  im¬ 
portant  property  of  this  approach  is  that  the  prediction  at  an  observation 
point  yields  the  observation  itself  (here  an  errorless  quantity)  with  the 
zero  prediction  error.  At  other  points  the  predictions  are  optimized  in  the 
sense  that  they  result  in  the  smallest  trace  of  the  dispersion  matrix  obtained 
using  the  M-averaging. 

In  practice,  the  covariance  function  may  be  computed  from  some  approxi¬ 
mate  values  of  the  S.H.  coefficients  complete  through  a  certain  degree  and 
order.  It  can  be  argued  that  this  function  should  be  compatible  with  the 
observations  and  should  ideally  be  computed  from  them,  since  it  seems  only 
reasonable  to  assume  that  if  one  is  making  predictions  which  are  based  on  some 
measured  values,  the  covariance  function  should  be  closely  related  to  them. 
However,  it  is  quite  possible  that  the  predicted  values  are  relatively  insensi¬ 
tive  to  variations  in  the  covariance  function.  Be  that  as  it  may,  caution 
should  be  exercised  when  interpreting  the  predicted  values  and  certain  predic- 
tion  errors  represented  by  Mice  }  .  The  covariance  function  as  used  in  these 
computations  expresses  a  general  behavior  of  the  theoretical  function  with 
respect  to  its  whole  domain,  the  earth.  Therefore,  individual  predictions  may 
be  somewhat  "blurred"  by  this  global  property,  while  being  "best"  overall. 

This  "blurring"  effect  is  even  more  noticeable  when  the  average  matrix  M{eeTl 
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is  considered,  since  e  is  expressed  in  terms  of  AT  which  already  depends  on 
the  covariance  function  and,  further,  since  the  product  ee  is  averaged  over 
the  whole  earth. 

An  important  practical  advantage  of  the  best  linear  prediction  method 
resides  in  avoiding  inversions  of  large  matrices.  For  the  price  of  introduc¬ 
ing  certain  (admissible)  errors,  only  a  small  number  of  observations  may  be 
used  to  predict  a  function  at  any  desired  location.  The  observatio*  ,’nts 

included  should  be  chosen  in  the  neighborhood  of  the  prediction  poi  The 
size  of  the  neighborhood  and  the  number  (n)  of  observation  points  ;"b,ject 
to  judgment.  The  matrix  H  to  be  inverted,  then,  has  the  dimensions  ,.i*n)  and 
most  often  is  likely  to  be  comfortably  small.  This  feature  is  attributed  to 
the  fact  that  the  covariance  functions  considered  decrease  rather  rapidly  with 
increasing  distances  from  prediction  points  since  geoid  undulations  and, 
especially,  gravity  anomalies  are  most  affected  by  local  and  regional  factors. 


APPENDIX  3 

TRANSLOCATION  ALGORITHM  ADAPTED  FOR  THE 
LEAST-SQUARES  COLLOCATION  WITH  NOISE 

In  reference  to  Section  6.4,  the  task  at  hand  can  be  presented  as 
follows.  Given  the  prediction  point  J  and  two  observation  points  numbered 
for  convenience  1  and  2,  find  the  seven  values  of  the  "storage"  matrix  for 
the  averaged  point  P  replacing  both  of  the  original  observation  points. 

These  points  are  considered  to  be  averaged  for  the  first  time  so  that  their 
weights  are  taken  as 

pi  =  ’  P2  =  ly^J2  *  ( A3,  la  ,b) 

If  they  were  computed,  the  coordinates  ( <J> » A )  of  P  would  be  given  by 

4>  =  (P14»1  +  P2$2)/(P1  +  P2)  »  X  *  (P1X1  +  P2^2)'(P,  +  P2)  1  (A3. 2a ,b) 

the  corresponding  value  of  the  geoid  undulation  at  P  is 

N  =  (piN1  +  p2N2)/(p1  +  p2)  •  (A3. 3) 

If  one  denotes 

d<f>  =  4»  -  ^  ,  dX  =  A2  -  ,  (A3.4a,b) 

it  follows  from  (A3.2a,b)  that 

(A3.5a,b,c) 

d<p  =  <t>  -  <1»1  =  w  d<|>  ,  dX1  =  A  -  \1  =  w  dX  ,  w  =  P2/(P1  +  P2)  • 
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In  fact,  upon  considering  a  (small)  rectangular  triangle  formed,  on  a  unit 
sphere,  by  the  segment  of  the  great  circle  between  1  and  2  and  by  the  cor¬ 
responding  segments  of  one  meridian  and  one  parallel,  one  would  also  deduce 


dd>1/d<p2  =  dX1/dX2  =  ty'/ty"  =  p2/p1  =  <j>J1/(|/J2  , 


where  d<|>2  =  4>  -  <J>  ,  dX2  =  X2  -  X  and  ip‘  ,  ip"  are  the  angular  distances  1-P 
and  P-2,  respectively,  measured  along  1-2.  In  other  words,  ip‘/<P 11  = 
gives  the  location  of  P  with  respect  to  the  locations  of  1  and  2  according  to 
the  desired  emphasis  on  the  distances  of  the  observation  points  1  and  2  from 
the  prediction  point  J. 

Next  designate  4>1  or  Xx  by  the  general  symbol  o^,  $2  or  X2  by  the  general 
symbol  a2,  and  (j>  or  X  by  the  general  symbol  a.  Consider  the  most  conservative 
case  when  dealing  with  round-off  errors  and  their  dependence  on  and  include 
second-order  effects  so  that  the  error  is  of  the  order  (o^-c^)3  at  the  most. 
Upon  applying  the  first  parts  of  (A3.5a,b)  to  "a",  one  has 


1  2 

sina  =  sino^  +  coso^do^  -  tj-  sina^do^  , 
cosa  =  coso^  -  sino^da^  -  ^  coso^da3  . 


(A3. 6a) 
(A3. 6b) 


The  following  relations,  which  should  be  understood  as  approximations  to  within 
the  desired  accuracy,  will  be  used: 

da  =  (sina2  -  sina^  /cos  a  ...  if  fcosaj  >  IsinaJ  ,  (A3. 7a) 

y  y  y 

da  =  (cosa2  -  cosa1)  /sin  a  ...  if  jcosaj  <  IsinaJ  .  (A3. 7b) 
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The  above  distinction  with  regard  to  a1  gives  rise  to  the  cases  "a"  and  "b" 
in  the  same  order. 

If  (A3.4a,b)  are  applied  to  "a",  one  obtains  (A3.6a,b)  with  a2  replacing 
a  and  da  replacing  dar  From  these  new  equations,  respectively,  and  from 
(A3.7a,b),  respectively,  it  then  follows  for  the  cases  a  and  b: 

da  =  T1(l  +  •^•tga1xT1)  ,  T1  =  (sina9  -  sina] )/cosaj  ,  (A3. 8a) 

da  =  -T^(l  +  ^cotga^T^)  ,  =  (cosa2  -  cosa^/sinaj  .  (A3. 8b) 

The  second  parts  of  (A3.5a,b)  applied  to  "a"  yield 

2  2  2 

da.^  =  w  da  ,  da1  =  w  da  , 

where  w  is  given  in  (A3. 5c),  da  is  given  in  (A3. 8a)  or  (A3. 8b),  and  da  is 
adopted  from  (A3. 7a)  or  (A3. 7b).  With  these  quantities,  equations  (A3.6a,b) 
are  readily  transcribed  as 

(A3. 9a) 

I  2  I 

sina  =  sinaj^  +  T2cosa1  -  ^T  sina1  ,  T  =wTJ  ,  T2  =  T  (l+^tga^  ; 

(A3. 9b) 

1  2  1 
cosa  =  cosaj^  +  T2sina1  -  jT'  cosa1  ,  T'  =  wTl[  ,  T2  =  T'  (1  +  -gTJcotga1)  . 

The  algorithm  can  now  be  summarized. 
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If  IcosaJ  >  |sina1|  : 


1  2 

sina  =  sina1  +  T2cosa1  -  -^-T  sina^^  ; 

2 

cosa  =  +(l-sin  a) 2  ,  sign(cosa)  =  signUosa.^  ; 

Tl  =  (sina2  -  sina1)/cosa1  ,  T  =  w  T1  ,  T2  =  T(  1  +  ^-T1sina1/cosa1)  . 

If  |cosa1|  <  IsinaJ  : 

1  2 

cosa  =  cosa1  +  T2sina1  -  -^T'  cosa1  ; 

2  V* 

sina  =  +(l-cos  a) 2  ,  sign(sina)  =  sign(sina1)  ; 

=  (cosa2  -  cosa1)/sina1  ,  T'=  wT|  ,  T2  =  T' (1 +|-T^cosa1/sina1)  . 

N  =  (p1N1  +  P2N2)/(Pl  +  P2)  • 

The  quantity  w  is  given  in  (A3. 5c)  while  p  and  p2  are  given  in  (A3.1a,b), 
respectively.  As  stipulated,  only  sine  and  cosine  functions  (needed  in  the 
"storage"  matrix  of  Section  6.4)  appear  in  this  algorithm.  With  regard  to  the 
angles  ,  these  functions  are  computed  by  the  standard  formulas. 
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